{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute health expenditures per income level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import stats \n",
    "from importlib import reload\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import os\n",
    "import sys\n",
    "from statsmodels.formula.api import ols,wls\n",
    "module_path = os.path.abspath(os.path.join('..'))\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from model import micro,macro, params\n",
    "from model import distributions as dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulated Gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "file = open('output/eq_pus_ge3.pkl','rb')\n",
    "eq = pickle.load(file)\n",
    "file.close()\n",
    "opt = pd.read_pickle('output/opt_pus_ge3.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wmean(x,var):\n",
    "    xx = x.loc[~x[var].isna(),:]\n",
    "    names = {var: (xx[var] * xx['ps']).sum()/xx['ps'].sum()} \n",
    "    return pd.Series(names, index=[var])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"10\" halign=\"left\">m</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>e</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>h</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.525576</td>\n",
       "      <td>7.810403e-01</td>\n",
       "      <td>1.29577</td>\n",
       "      <td>2.689467</td>\n",
       "      <td>7.29176</td>\n",
       "      <td>13.616625</td>\n",
       "      <td>16.92276</td>\n",
       "      <td>18.823232</td>\n",
       "      <td>2.058494e+01</td>\n",
       "      <td>22.230655</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>5.656805e-18</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>4.268590e-14</td>\n",
       "      <td>0.019326</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          m                                                                 \\\n",
       "e         0             1        2         3        4          5         6   \n",
       "h                                                                            \n",
       "0  0.525576  7.810403e-01  1.29577  2.689467  7.29176  13.616625  16.92276   \n",
       "1  0.000000  5.656805e-18  0.00000  0.000000  0.00000   0.000000   0.00000   \n",
       "\n",
       "                                       \n",
       "e          7             8          9  \n",
       "h                                      \n",
       "0  18.823232  2.058494e+01  22.230655  \n",
       "1   0.000000  4.268590e-14   0.019326  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tab_eh = opt.groupby(['h','e']).apply(wmean,var='m').unstack()\n",
    "tab_eh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   e\n",
       "h  0    0.677376\n",
       "   1    0.719265\n",
       "   2    0.778303\n",
       "   3    0.846530\n",
       "   4    0.910622\n",
       "   5    0.949463\n",
       "   6    0.962733\n",
       "   7    0.965976\n",
       "   8    0.967188\n",
       "   9    0.967949\n",
       "dtype: float64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tab_h = opt.groupby('e').apply(wmean,var='h').unstack()\n",
    "tab_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>m</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>e</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.169564</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.219266</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.287269</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.412752</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.651721</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.688138</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.630657</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.640445</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.675423</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.731230</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          m\n",
       "e          \n",
       "0  0.169564\n",
       "1  0.219266\n",
       "2  0.287269\n",
       "3  0.412752\n",
       "4  0.651721\n",
       "5  0.688138\n",
       "6  0.630657\n",
       "7  0.640445\n",
       "8  0.675423\n",
       "9  0.731230"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_m_sim = opt.groupby('e').apply(wmean,var='m')\n",
    "grad_m_sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>m</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-0.739822</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.663559</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.559215</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-0.366673</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.055878</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>-0.032321</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>-0.017302</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.036369</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.121999</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           m\n",
       "1  -0.739822\n",
       "2  -0.663559\n",
       "3  -0.559215\n",
       "4  -0.366673\n",
       "5   0.000000\n",
       "6   0.055878\n",
       "7  -0.032321\n",
       "8  -0.017302\n",
       "9   0.036369\n",
       "10  0.121999"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_m_sim['m'] = grad_m_sim['m']/grad_m_sim.loc[4,'m']-1\n",
    "grad_m_sim.index = [x for x in range(1,11)]\n",
    "grad_m_sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>duid</th>\n",
       "      <th>dupersid</th>\n",
       "      <th>famidyr</th>\n",
       "      <th>famszeyr</th>\n",
       "      <th>famrfpyr</th>\n",
       "      <th>region53</th>\n",
       "      <th>age</th>\n",
       "      <th>racex</th>\n",
       "      <th>hispan</th>\n",
       "      <th>marry53x</th>\n",
       "      <th>...</th>\n",
       "      <th>nh_diabe</th>\n",
       "      <th>nh_lunge</th>\n",
       "      <th>nh_cancre</th>\n",
       "      <th>smokev</th>\n",
       "      <th>nh_smoken</th>\n",
       "      <th>mergeMEPSNHIS</th>\n",
       "      <th>smoken</th>\n",
       "      <th>totinc</th>\n",
       "      <th>inscov</th>\n",
       "      <th>_merge</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>20002.0</td>\n",
       "      <td>20002014</td>\n",
       "      <td>A</td>\n",
       "      <td>6.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>30000.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>20002.0</td>\n",
       "      <td>20002014</td>\n",
       "      <td>A</td>\n",
       "      <td>6.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>25200.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>20002.0</td>\n",
       "      <td>20002021</td>\n",
       "      <td>A</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>20002.0</td>\n",
       "      <td>20002021</td>\n",
       "      <td>A</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>20002.0</td>\n",
       "      <td>20002038</td>\n",
       "      <td>A</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 125 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      duid  dupersid famidyr  famszeyr  famrfpyr  region53   age  racex  \\\n",
       "0  20002.0  20002014       A       6.0       1.0       4.0  39.0    1.0   \n",
       "1  20002.0  20002014       A       6.0       1.0       4.0  40.0    1.0   \n",
       "2  20002.0  20002021       A       6.0       0.0       4.0  16.0    1.0   \n",
       "3  20002.0  20002021       A       6.0       0.0       4.0  17.0    1.0   \n",
       "4  20002.0  20002038       A       6.0       0.0       4.0  14.0    1.0   \n",
       "\n",
       "   hispan  marry53x  ...  nh_diabe  nh_lunge  nh_cancre  smokev  nh_smoken  \\\n",
       "0     1.0       1.0  ...       NaN       NaN        NaN     NaN        NaN   \n",
       "1     1.0       1.0  ...       NaN       NaN        NaN     NaN        NaN   \n",
       "2     1.0       5.0  ...       NaN       NaN        NaN     NaN        NaN   \n",
       "3     1.0       5.0  ...       NaN       NaN        NaN     NaN        NaN   \n",
       "4     1.0       6.0  ...       NaN       NaN        NaN     NaN        NaN   \n",
       "\n",
       "   mergeMEPSNHIS  smoken   totinc  inscov  _merge  \n",
       "0            1.0     0.0  30000.0     1.0       3  \n",
       "1            1.0     0.0  25200.0     1.0       3  \n",
       "2            1.0     NaN      0.0     3.0       3  \n",
       "3            1.0     NaN      0.0     3.0       3  \n",
       "4            1.0     NaN      0.0     3.0       3  \n",
       "\n",
       "[5 rows x 125 columns]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_stata('../data_sources/meps/MEPS0004_agghealth.dta', convert_categoricals=False)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['duid',\n",
       " 'dupersid',\n",
       " 'famidyr',\n",
       " 'famszeyr',\n",
       " 'famrfpyr',\n",
       " 'region53',\n",
       " 'age',\n",
       " 'racex',\n",
       " 'hispan',\n",
       " 'marry53x',\n",
       " 'educyear',\n",
       " 'totded04',\n",
       " 'rthlth53',\n",
       " 'mnhlth53',\n",
       " 'iadlhp53',\n",
       " 'adlhlp53',\n",
       " 'bmindx53',\n",
       " 'diabdx53',\n",
       " 'asthdx53',\n",
       " 'hibpdx53',\n",
       " 'chddx53',\n",
       " 'angidx53',\n",
       " 'midx53',\n",
       " 'ohrtdx53',\n",
       " 'strkdx53',\n",
       " 'emphdx53',\n",
       " 'adsmok42',\n",
       " 'tottch04',\n",
       " 'totexp',\n",
       " 'totslf',\n",
       " 'totmcr',\n",
       " 'totmcd',\n",
       " 'totprv',\n",
       " 'totva',\n",
       " 'tottri04',\n",
       " 'totofd',\n",
       " 'totstl04',\n",
       " 'totwcp',\n",
       " 'totopr',\n",
       " 'totopu',\n",
       " 'totosr',\n",
       " 'doctim',\n",
       " 'hsptim',\n",
       " 'hspnit',\n",
       " 'perwt',\n",
       " 'famwt',\n",
       " 'male',\n",
       " 'mepsexp',\n",
       " 'mepsslf',\n",
       " 'mepsmcr',\n",
       " 'mepsmcd',\n",
       " 'mepsprv',\n",
       " 'mepsva',\n",
       " 'mepsofd',\n",
       " 'mepswcp',\n",
       " 'mepsopr',\n",
       " 'mepsopu',\n",
       " 'mepsosr',\n",
       " 'yr',\n",
       " 'totded03',\n",
       " 'tottch03',\n",
       " 'tottri03',\n",
       " 'totstl03',\n",
       " 'totded02',\n",
       " 'tottch02',\n",
       " 'tottri02',\n",
       " 'totstl02',\n",
       " 'totded01',\n",
       " 'tottch01',\n",
       " 'tottri01',\n",
       " 'totstl01',\n",
       " 'totded00',\n",
       " 'tottch00',\n",
       " 'tottri00',\n",
       " 'totstl00',\n",
       " 'cccodex',\n",
       " 'ncccodex',\n",
       " 'diabecr',\n",
       " 'hibpecr',\n",
       " 'strokecr',\n",
       " 'heartecr',\n",
       " 'lungecr',\n",
       " 'cancrecr',\n",
       " 'hearte',\n",
       " 'diabe',\n",
       " 'hibpe',\n",
       " 'stroke',\n",
       " 'lunge',\n",
       " 'adl1p',\n",
       " 'iadl1',\n",
       " 'age5659',\n",
       " 'age6064',\n",
       " 'black',\n",
       " 'educ',\n",
       " 'hsdrop',\n",
       " 'somecol',\n",
       " 'collgrad',\n",
       " 'regnth',\n",
       " 'regmid',\n",
       " 'regwst',\n",
       " 'widowed',\n",
       " 'single',\n",
       " 'bmi',\n",
       " 'obese',\n",
       " 'overwt',\n",
       " 'normalwt',\n",
       " 'underwt',\n",
       " 'wtstate',\n",
       " 'nh_hhx',\n",
       " 'nh_px',\n",
       " 'panel',\n",
       " 'nhisyr',\n",
       " 'nh_hibpe',\n",
       " 'nh_hearte',\n",
       " 'nh_stroke',\n",
       " 'nh_diabe',\n",
       " 'nh_lunge',\n",
       " 'nh_cancre',\n",
       " 'smokev',\n",
       " 'nh_smoken',\n",
       " 'mergeMEPSNHIS',\n",
       " 'smoken',\n",
       " 'totinc',\n",
       " 'inscov',\n",
       " '_merge']"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.columns.to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[df['yr'].isin([2003,2004,2005,2006,2007])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[(df.age>=35)&(df.age<=75)]\n",
    "df = df[df.totinc>10e3]\n",
    "df = df[~df.mepsexp.isna()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = df[['totinc','adl1p','age','mepsexp','totmcd','yr','male','obese','diabe','hearte','hibpe','stroke','lunge','perwt','smoken','smokev']].dropna(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['loginc'] = np.log(data['totinc'])\n",
    "data['m'] = data['mepsexp'] \n",
    "data['logm'] = np.log(data.loc[data.m>0,'m'].min() + data['m'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "42.62466796875"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.01*data['m'].mean() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/loulou/anaconda3/lib/python3.8/site-packages/statsmodels/base/model.py:127: ValueWarning: unknown kwargs ['var_weights']\n",
      "  warnings.warn(msg, ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "model = wls('loginc ~ C(age) + C(yr) ',data=data,var_weights=np.asarray(data['perwt'])).fit()\n",
    "# + smoken + smokev + obese + diabe + stroke + hibpe + lunge + hearte + male \n",
    "data.loc[:,'ey'] = model.resid.to_list()\n",
    "data.ey = (data.ey - data.ey.mean())/data.ey.std()\n",
    "data['ey'] = data['ey'].clip(lower=-2.5,upper=2.5)\n",
    "data['ey10'] = pd.cut(data['ey'],bins=10,labels=[x for x in range(1,11)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "for q in range(1,11):\n",
    "\tdata['ey10_'+str(q)] = np.where(data['ey10']==q,1,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>totinc</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>39530.646980</td>\n",
       "      <td>29307.996623</td>\n",
       "      <td>10001.000000</td>\n",
       "      <td>19410.500000</td>\n",
       "      <td>30834.500000</td>\n",
       "      <td>50118.750000</td>\n",
       "      <td>279668.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>adl1p</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.011678</td>\n",
       "      <td>0.107438</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>51.375365</td>\n",
       "      <td>11.167422</td>\n",
       "      <td>35.000000</td>\n",
       "      <td>42.000000</td>\n",
       "      <td>50.000000</td>\n",
       "      <td>60.000000</td>\n",
       "      <td>75.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mepsexp</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>4262.466797</td>\n",
       "      <td>9039.359375</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>361.830154</td>\n",
       "      <td>1484.150269</td>\n",
       "      <td>4296.375122</td>\n",
       "      <td>205658.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>totmcd</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>260.089939</td>\n",
       "      <td>2611.798043</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>113045.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>yr</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>2003.534302</td>\n",
       "      <td>0.500991</td>\n",
       "      <td>2003.000000</td>\n",
       "      <td>2003.000000</td>\n",
       "      <td>2004.000000</td>\n",
       "      <td>2004.000000</td>\n",
       "      <td>2004.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>male</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.458567</td>\n",
       "      <td>0.498305</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>obese</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.299349</td>\n",
       "      <td>0.457999</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diabe</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.084101</td>\n",
       "      <td>0.277563</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hearte</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.104424</td>\n",
       "      <td>0.305836</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hibpe</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.310128</td>\n",
       "      <td>0.462586</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>stroke</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.024590</td>\n",
       "      <td>0.154888</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>lunge</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.017965</td>\n",
       "      <td>0.132829</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>perwt</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>10825.980579</td>\n",
       "      <td>5811.769553</td>\n",
       "      <td>711.079529</td>\n",
       "      <td>6759.226783</td>\n",
       "      <td>10100.714239</td>\n",
       "      <td>13874.193993</td>\n",
       "      <td>57180.156250</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>smoken</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.222883</td>\n",
       "      <td>0.416204</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>smokev</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.482596</td>\n",
       "      <td>0.499732</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>loginc</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>10.368490</td>\n",
       "      <td>0.642292</td>\n",
       "      <td>9.210440</td>\n",
       "      <td>9.873569</td>\n",
       "      <td>10.336389</td>\n",
       "      <td>10.822150</td>\n",
       "      <td>12.541358</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>m</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>4262.466797</td>\n",
       "      <td>9039.359375</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>361.830154</td>\n",
       "      <td>1484.150269</td>\n",
       "      <td>4296.375122</td>\n",
       "      <td>205658.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>logm</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>6.714816</td>\n",
       "      <td>2.538181</td>\n",
       "      <td>0.734520</td>\n",
       "      <td>5.896919</td>\n",
       "      <td>7.304001</td>\n",
       "      <td>8.366012</td>\n",
       "      <td>12.233980</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>-0.001244</td>\n",
       "      <td>0.996690</td>\n",
       "      <td>-2.028125</td>\n",
       "      <td>-0.764337</td>\n",
       "      <td>-0.054477</td>\n",
       "      <td>0.699739</td>\n",
       "      <td>2.500000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_1</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.045587</td>\n",
       "      <td>0.208600</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_2</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.098585</td>\n",
       "      <td>0.298121</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_3</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.137660</td>\n",
       "      <td>0.344562</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_4</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.159780</td>\n",
       "      <td>0.366422</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_5</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.163822</td>\n",
       "      <td>0.370135</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_6</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.141815</td>\n",
       "      <td>0.348879</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_7</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.117224</td>\n",
       "      <td>0.321705</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_8</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.075118</td>\n",
       "      <td>0.263596</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_9</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.028520</td>\n",
       "      <td>0.166463</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ey10_10</th>\n",
       "      <td>8906.0</td>\n",
       "      <td>0.031889</td>\n",
       "      <td>0.175713</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          count          mean           std           min           25%  \\\n",
       "totinc   8906.0  39530.646980  29307.996623  10001.000000  19410.500000   \n",
       "adl1p    8906.0      0.011678      0.107438      0.000000      0.000000   \n",
       "age      8906.0     51.375365     11.167422     35.000000     42.000000   \n",
       "mepsexp  8906.0   4262.466797   9039.359375      0.000000    361.830154   \n",
       "totmcd   8906.0    260.089939   2611.798043      0.000000      0.000000   \n",
       "yr       8906.0   2003.534302      0.500991   2003.000000   2003.000000   \n",
       "male     8906.0      0.458567      0.498305      0.000000      0.000000   \n",
       "obese    8906.0      0.299349      0.457999      0.000000      0.000000   \n",
       "diabe    8906.0      0.084101      0.277563      0.000000      0.000000   \n",
       "hearte   8906.0      0.104424      0.305836      0.000000      0.000000   \n",
       "hibpe    8906.0      0.310128      0.462586      0.000000      0.000000   \n",
       "stroke   8906.0      0.024590      0.154888      0.000000      0.000000   \n",
       "lunge    8906.0      0.017965      0.132829      0.000000      0.000000   \n",
       "perwt    8906.0  10825.980579   5811.769553    711.079529   6759.226783   \n",
       "smoken   8906.0      0.222883      0.416204      0.000000      0.000000   \n",
       "smokev   8906.0      0.482596      0.499732      0.000000      0.000000   \n",
       "loginc   8906.0     10.368490      0.642292      9.210440      9.873569   \n",
       "m        8906.0   4262.466797   9039.359375      0.000000    361.830154   \n",
       "logm     8906.0      6.714816      2.538181      0.734520      5.896919   \n",
       "ey       8906.0     -0.001244      0.996690     -2.028125     -0.764337   \n",
       "ey10_1   8906.0      0.045587      0.208600      0.000000      0.000000   \n",
       "ey10_2   8906.0      0.098585      0.298121      0.000000      0.000000   \n",
       "ey10_3   8906.0      0.137660      0.344562      0.000000      0.000000   \n",
       "ey10_4   8906.0      0.159780      0.366422      0.000000      0.000000   \n",
       "ey10_5   8906.0      0.163822      0.370135      0.000000      0.000000   \n",
       "ey10_6   8906.0      0.141815      0.348879      0.000000      0.000000   \n",
       "ey10_7   8906.0      0.117224      0.321705      0.000000      0.000000   \n",
       "ey10_8   8906.0      0.075118      0.263596      0.000000      0.000000   \n",
       "ey10_9   8906.0      0.028520      0.166463      0.000000      0.000000   \n",
       "ey10_10  8906.0      0.031889      0.175713      0.000000      0.000000   \n",
       "\n",
       "                  50%           75%            max  \n",
       "totinc   30834.500000  50118.750000  279668.000000  \n",
       "adl1p        0.000000      0.000000       1.000000  \n",
       "age         50.000000     60.000000      75.000000  \n",
       "mepsexp   1484.150269   4296.375122  205658.000000  \n",
       "totmcd       0.000000      0.000000  113045.000000  \n",
       "yr        2004.000000   2004.000000    2004.000000  \n",
       "male         0.000000      1.000000       1.000000  \n",
       "obese        0.000000      1.000000       1.000000  \n",
       "diabe        0.000000      0.000000       1.000000  \n",
       "hearte       0.000000      0.000000       1.000000  \n",
       "hibpe        0.000000      1.000000       1.000000  \n",
       "stroke       0.000000      0.000000       1.000000  \n",
       "lunge        0.000000      0.000000       1.000000  \n",
       "perwt    10100.714239  13874.193993   57180.156250  \n",
       "smoken       0.000000      0.000000       1.000000  \n",
       "smokev       0.000000      1.000000       1.000000  \n",
       "loginc      10.336389     10.822150      12.541358  \n",
       "m         1484.150269   4296.375122  205658.000000  \n",
       "logm         7.304001      8.366012      12.233980  \n",
       "ey          -0.054477      0.699739       2.500000  \n",
       "ey10_1       0.000000      0.000000       1.000000  \n",
       "ey10_2       0.000000      0.000000       1.000000  \n",
       "ey10_3       0.000000      0.000000       1.000000  \n",
       "ey10_4       0.000000      0.000000       1.000000  \n",
       "ey10_5       0.000000      0.000000       1.000000  \n",
       "ey10_6       0.000000      0.000000       1.000000  \n",
       "ey10_7       0.000000      0.000000       1.000000  \n",
       "ey10_8       0.000000      0.000000       1.000000  \n",
       "ey10_9       0.000000      0.000000       1.000000  \n",
       "ey10_10      0.000000      0.000000       1.000000  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.describe().transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/loulou/anaconda3/lib/python3.8/site-packages/statsmodels/base/model.py:127: ValueWarning: unknown kwargs ['var_weights']\n",
      "  warnings.warn(msg, ValueWarning)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>WLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>logm</td>       <th>  R-squared:         </th> <td>   0.113</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>WLS</td>       <th>  Adj. R-squared:    </th> <td>   0.108</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   22.53</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 17 May 2022</td> <th>  Prob (F-statistic):</th> <td>5.21e-190</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>11:14:04</td>     <th>  Log-Likelihood:    </th> <td> -20399.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  8906</td>      <th>  AIC:               </th> <td>4.090e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  8855</td>      <th>  BIC:               </th> <td>4.126e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>    50</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>    5.7718</td> <td>    0.158</td> <td>   36.427</td> <td> 0.000</td> <td>    5.461</td> <td>    6.082</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.36.0]</th>  <td>    0.0733</td> <td>    0.200</td> <td>    0.366</td> <td> 0.714</td> <td>   -0.319</td> <td>    0.466</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.37.0]</th>  <td>   -0.1080</td> <td>    0.202</td> <td>   -0.535</td> <td> 0.593</td> <td>   -0.504</td> <td>    0.288</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.38.0]</th>  <td>    0.0320</td> <td>    0.201</td> <td>    0.159</td> <td> 0.873</td> <td>   -0.362</td> <td>    0.426</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.39.0]</th>  <td>   -0.0883</td> <td>    0.201</td> <td>   -0.439</td> <td> 0.661</td> <td>   -0.483</td> <td>    0.306</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.40.0]</th>  <td>    0.2469</td> <td>    0.199</td> <td>    1.239</td> <td> 0.215</td> <td>   -0.144</td> <td>    0.637</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.41.0]</th>  <td>    0.2374</td> <td>    0.195</td> <td>    1.215</td> <td> 0.224</td> <td>   -0.146</td> <td>    0.620</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.42.0]</th>  <td>    0.2928</td> <td>    0.203</td> <td>    1.444</td> <td> 0.149</td> <td>   -0.105</td> <td>    0.690</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.43.0]</th>  <td>    0.2471</td> <td>    0.204</td> <td>    1.209</td> <td> 0.227</td> <td>   -0.153</td> <td>    0.648</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.44.0]</th>  <td>    0.5108</td> <td>    0.203</td> <td>    2.520</td> <td> 0.012</td> <td>    0.114</td> <td>    0.908</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.45.0]</th>  <td>    0.8884</td> <td>    0.200</td> <td>    4.450</td> <td> 0.000</td> <td>    0.497</td> <td>    1.280</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.46.0]</th>  <td>    0.8845</td> <td>    0.208</td> <td>    4.244</td> <td> 0.000</td> <td>    0.476</td> <td>    1.293</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.47.0]</th>  <td>    0.6933</td> <td>    0.203</td> <td>    3.414</td> <td> 0.001</td> <td>    0.295</td> <td>    1.091</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.48.0]</th>  <td>    0.7154</td> <td>    0.208</td> <td>    3.444</td> <td> 0.001</td> <td>    0.308</td> <td>    1.123</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.49.0]</th>  <td>    0.8715</td> <td>    0.208</td> <td>    4.191</td> <td> 0.000</td> <td>    0.464</td> <td>    1.279</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.50.0]</th>  <td>    0.9586</td> <td>    0.207</td> <td>    4.638</td> <td> 0.000</td> <td>    0.553</td> <td>    1.364</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.51.0]</th>  <td>    1.1212</td> <td>    0.211</td> <td>    5.302</td> <td> 0.000</td> <td>    0.707</td> <td>    1.536</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.52.0]</th>  <td>    1.1587</td> <td>    0.208</td> <td>    5.576</td> <td> 0.000</td> <td>    0.751</td> <td>    1.566</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.53.0]</th>  <td>    1.1226</td> <td>    0.215</td> <td>    5.211</td> <td> 0.000</td> <td>    0.700</td> <td>    1.545</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.54.0]</th>  <td>    1.3132</td> <td>    0.213</td> <td>    6.154</td> <td> 0.000</td> <td>    0.895</td> <td>    1.731</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.55.0]</th>  <td>    1.4092</td> <td>    0.217</td> <td>    6.491</td> <td> 0.000</td> <td>    0.984</td> <td>    1.835</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.56.0]</th>  <td>    1.1776</td> <td>    0.214</td> <td>    5.509</td> <td> 0.000</td> <td>    0.759</td> <td>    1.597</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.57.0]</th>  <td>    1.5739</td> <td>    0.218</td> <td>    7.232</td> <td> 0.000</td> <td>    1.147</td> <td>    2.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.58.0]</th>  <td>    1.4938</td> <td>    0.227</td> <td>    6.578</td> <td> 0.000</td> <td>    1.049</td> <td>    1.939</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.59.0]</th>  <td>    1.9002</td> <td>    0.239</td> <td>    7.956</td> <td> 0.000</td> <td>    1.432</td> <td>    2.368</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.60.0]</th>  <td>    1.6420</td> <td>    0.235</td> <td>    6.994</td> <td> 0.000</td> <td>    1.182</td> <td>    2.102</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.61.0]</th>  <td>    1.8924</td> <td>    0.236</td> <td>    8.011</td> <td> 0.000</td> <td>    1.429</td> <td>    2.355</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.62.0]</th>  <td>    2.0099</td> <td>    0.233</td> <td>    8.636</td> <td> 0.000</td> <td>    1.554</td> <td>    2.466</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.63.0]</th>  <td>    1.8638</td> <td>    0.243</td> <td>    7.661</td> <td> 0.000</td> <td>    1.387</td> <td>    2.341</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.64.0]</th>  <td>    1.8437</td> <td>    0.247</td> <td>    7.467</td> <td> 0.000</td> <td>    1.360</td> <td>    2.328</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.65.0]</th>  <td>    2.3343</td> <td>    0.248</td> <td>    9.413</td> <td> 0.000</td> <td>    1.848</td> <td>    2.820</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.66.0]</th>  <td>    1.9793</td> <td>    0.246</td> <td>    8.050</td> <td> 0.000</td> <td>    1.497</td> <td>    2.461</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.67.0]</th>  <td>    2.0425</td> <td>    0.257</td> <td>    7.935</td> <td> 0.000</td> <td>    1.538</td> <td>    2.547</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.68.0]</th>  <td>    2.2760</td> <td>    0.248</td> <td>    9.195</td> <td> 0.000</td> <td>    1.791</td> <td>    2.761</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.69.0]</th>  <td>    2.2431</td> <td>    0.259</td> <td>    8.644</td> <td> 0.000</td> <td>    1.734</td> <td>    2.752</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.70.0]</th>  <td>    2.5740</td> <td>    0.262</td> <td>    9.836</td> <td> 0.000</td> <td>    2.061</td> <td>    3.087</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.71.0]</th>  <td>    2.0607</td> <td>    0.268</td> <td>    7.692</td> <td> 0.000</td> <td>    1.536</td> <td>    2.586</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.72.0]</th>  <td>    2.6001</td> <td>    0.249</td> <td>   10.429</td> <td> 0.000</td> <td>    2.111</td> <td>    3.089</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.73.0]</th>  <td>    2.3910</td> <td>    0.250</td> <td>    9.567</td> <td> 0.000</td> <td>    1.901</td> <td>    2.881</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.74.0]</th>  <td>    2.6128</td> <td>    0.267</td> <td>    9.777</td> <td> 0.000</td> <td>    2.089</td> <td>    3.137</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.75.0]</th>  <td>    2.6292</td> <td>    0.258</td> <td>   10.182</td> <td> 0.000</td> <td>    2.123</td> <td>    3.135</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(yr)[T.2004.0]</th> <td>   -0.0472</td> <td>    0.051</td> <td>   -0.927</td> <td> 0.354</td> <td>   -0.147</td> <td>    0.053</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_1</th>          <td>   -0.3109</td> <td>    0.135</td> <td>   -2.295</td> <td> 0.022</td> <td>   -0.576</td> <td>   -0.045</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_2</th>          <td>   -0.5245</td> <td>    0.103</td> <td>   -5.110</td> <td> 0.000</td> <td>   -0.726</td> <td>   -0.323</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_3</th>          <td>   -0.4295</td> <td>    0.093</td> <td>   -4.609</td> <td> 0.000</td> <td>   -0.612</td> <td>   -0.247</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_4</th>          <td>   -0.2186</td> <td>    0.090</td> <td>   -2.440</td> <td> 0.015</td> <td>   -0.394</td> <td>   -0.043</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_6</th>          <td>    0.1822</td> <td>    0.092</td> <td>    1.973</td> <td> 0.048</td> <td>    0.001</td> <td>    0.363</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_7</th>          <td>    0.1934</td> <td>    0.097</td> <td>    1.986</td> <td> 0.047</td> <td>    0.002</td> <td>    0.384</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_8</th>          <td>    0.2131</td> <td>    0.112</td> <td>    1.900</td> <td> 0.057</td> <td>   -0.007</td> <td>    0.433</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_9</th>          <td>    0.1363</td> <td>    0.163</td> <td>    0.834</td> <td> 0.404</td> <td>   -0.184</td> <td>    0.456</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_10</th>         <td>    0.1874</td> <td>    0.156</td> <td>    1.202</td> <td> 0.229</td> <td>   -0.118</td> <td>    0.493</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>1307.653</td> <th>  Durbin-Watson:     </th> <td>   1.767</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>  <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1979.866</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>           <td>-1.065</td>  <th>  Prob(JB):          </th> <td>    0.00</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>       <td> 3.892</td>  <th>  Cond. No.          </th> <td>    44.8</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            WLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                   logm   R-squared:                       0.113\n",
       "Model:                            WLS   Adj. R-squared:                  0.108\n",
       "Method:                 Least Squares   F-statistic:                     22.53\n",
       "Date:                Tue, 17 May 2022   Prob (F-statistic):          5.21e-190\n",
       "Time:                        11:14:04   Log-Likelihood:                -20399.\n",
       "No. Observations:                8906   AIC:                         4.090e+04\n",
       "Df Residuals:                    8855   BIC:                         4.126e+04\n",
       "Df Model:                          50                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "===================================================================================\n",
       "                      coef    std err          t      P>|t|      [0.025      0.975]\n",
       "-----------------------------------------------------------------------------------\n",
       "Intercept           5.7718      0.158     36.427      0.000       5.461       6.082\n",
       "C(age)[T.36.0]      0.0733      0.200      0.366      0.714      -0.319       0.466\n",
       "C(age)[T.37.0]     -0.1080      0.202     -0.535      0.593      -0.504       0.288\n",
       "C(age)[T.38.0]      0.0320      0.201      0.159      0.873      -0.362       0.426\n",
       "C(age)[T.39.0]     -0.0883      0.201     -0.439      0.661      -0.483       0.306\n",
       "C(age)[T.40.0]      0.2469      0.199      1.239      0.215      -0.144       0.637\n",
       "C(age)[T.41.0]      0.2374      0.195      1.215      0.224      -0.146       0.620\n",
       "C(age)[T.42.0]      0.2928      0.203      1.444      0.149      -0.105       0.690\n",
       "C(age)[T.43.0]      0.2471      0.204      1.209      0.227      -0.153       0.648\n",
       "C(age)[T.44.0]      0.5108      0.203      2.520      0.012       0.114       0.908\n",
       "C(age)[T.45.0]      0.8884      0.200      4.450      0.000       0.497       1.280\n",
       "C(age)[T.46.0]      0.8845      0.208      4.244      0.000       0.476       1.293\n",
       "C(age)[T.47.0]      0.6933      0.203      3.414      0.001       0.295       1.091\n",
       "C(age)[T.48.0]      0.7154      0.208      3.444      0.001       0.308       1.123\n",
       "C(age)[T.49.0]      0.8715      0.208      4.191      0.000       0.464       1.279\n",
       "C(age)[T.50.0]      0.9586      0.207      4.638      0.000       0.553       1.364\n",
       "C(age)[T.51.0]      1.1212      0.211      5.302      0.000       0.707       1.536\n",
       "C(age)[T.52.0]      1.1587      0.208      5.576      0.000       0.751       1.566\n",
       "C(age)[T.53.0]      1.1226      0.215      5.211      0.000       0.700       1.545\n",
       "C(age)[T.54.0]      1.3132      0.213      6.154      0.000       0.895       1.731\n",
       "C(age)[T.55.0]      1.4092      0.217      6.491      0.000       0.984       1.835\n",
       "C(age)[T.56.0]      1.1776      0.214      5.509      0.000       0.759       1.597\n",
       "C(age)[T.57.0]      1.5739      0.218      7.232      0.000       1.147       2.000\n",
       "C(age)[T.58.0]      1.4938      0.227      6.578      0.000       1.049       1.939\n",
       "C(age)[T.59.0]      1.9002      0.239      7.956      0.000       1.432       2.368\n",
       "C(age)[T.60.0]      1.6420      0.235      6.994      0.000       1.182       2.102\n",
       "C(age)[T.61.0]      1.8924      0.236      8.011      0.000       1.429       2.355\n",
       "C(age)[T.62.0]      2.0099      0.233      8.636      0.000       1.554       2.466\n",
       "C(age)[T.63.0]      1.8638      0.243      7.661      0.000       1.387       2.341\n",
       "C(age)[T.64.0]      1.8437      0.247      7.467      0.000       1.360       2.328\n",
       "C(age)[T.65.0]      2.3343      0.248      9.413      0.000       1.848       2.820\n",
       "C(age)[T.66.0]      1.9793      0.246      8.050      0.000       1.497       2.461\n",
       "C(age)[T.67.0]      2.0425      0.257      7.935      0.000       1.538       2.547\n",
       "C(age)[T.68.0]      2.2760      0.248      9.195      0.000       1.791       2.761\n",
       "C(age)[T.69.0]      2.2431      0.259      8.644      0.000       1.734       2.752\n",
       "C(age)[T.70.0]      2.5740      0.262      9.836      0.000       2.061       3.087\n",
       "C(age)[T.71.0]      2.0607      0.268      7.692      0.000       1.536       2.586\n",
       "C(age)[T.72.0]      2.6001      0.249     10.429      0.000       2.111       3.089\n",
       "C(age)[T.73.0]      2.3910      0.250      9.567      0.000       1.901       2.881\n",
       "C(age)[T.74.0]      2.6128      0.267      9.777      0.000       2.089       3.137\n",
       "C(age)[T.75.0]      2.6292      0.258     10.182      0.000       2.123       3.135\n",
       "C(yr)[T.2004.0]    -0.0472      0.051     -0.927      0.354      -0.147       0.053\n",
       "ey10_1             -0.3109      0.135     -2.295      0.022      -0.576      -0.045\n",
       "ey10_2             -0.5245      0.103     -5.110      0.000      -0.726      -0.323\n",
       "ey10_3             -0.4295      0.093     -4.609      0.000      -0.612      -0.247\n",
       "ey10_4             -0.2186      0.090     -2.440      0.015      -0.394      -0.043\n",
       "ey10_6              0.1822      0.092      1.973      0.048       0.001       0.363\n",
       "ey10_7              0.1934      0.097      1.986      0.047       0.002       0.384\n",
       "ey10_8              0.2131      0.112      1.900      0.057      -0.007       0.433\n",
       "ey10_9              0.1363      0.163      0.834      0.404      -0.184       0.456\n",
       "ey10_10             0.1874      0.156      1.202      0.229      -0.118       0.493\n",
       "==============================================================================\n",
       "Omnibus:                     1307.653   Durbin-Watson:                   1.767\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1979.866\n",
       "Skew:                          -1.065   Prob(JB):                         0.00\n",
       "Kurtosis:                       3.892   Cond. No.                         44.8\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = wls('logm ~ ey10_1 +ey10_2 + ey10_3 + ey10_4 + ey10_6 + ey10_7 + ey10_8 + ey10_9 + ey10_10 + C(age) + C(yr) ',data=data,var_weights=np.asarray(data['perwt'])).fit()\n",
    "#+ smoken + smokev + obese + diabe + stroke + hibpe + lunge + hearte + male\n",
    "model.summary()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1     0.156416\n",
       "2     0.114373\n",
       "3     0.095196\n",
       "4     0.090118\n",
       "5     0.000000\n",
       "6     0.085721\n",
       "7     0.090745\n",
       "8     0.100431\n",
       "9     0.127536\n",
       "10    0.131535\n",
       "dtype: float64"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_m_data_se = model.HC0_se[['ey10_'+str(q) for q in range(1,11) if q!=5]]\n",
    "grad_m_data_se.index = [x for x in range(1,11) if x!=5]\n",
    "grad_m_data_se[5] = 0\n",
    "grad_m_data_se = grad_m_data_se.sort_index()\n",
    "grad_m_data_se"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    -0.310938\n",
       "2    -0.524490\n",
       "3    -0.429465\n",
       "4    -0.218555\n",
       "5     0.000000\n",
       "6     0.182247\n",
       "7     0.193383\n",
       "8     0.213062\n",
       "9     0.136299\n",
       "10    0.187443\n",
       "dtype: float64"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_m_data = model.params[['ey10_'+str(q) for q in range(1,11) if q!=5]]\n",
    "grad_m_data.index = [x for x in range(1,11) if x!=5]\n",
    "grad_m_data[5] = 0\n",
    "grad_m_data = grad_m_data.sort_index()\n",
    "grad_m_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/loulou/anaconda3/lib/python3.8/site-packages/statsmodels/base/model.py:127: ValueWarning: unknown kwargs ['var_weights']\n",
      "  warnings.warn(msg, ValueWarning)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>WLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>logm</td>       <th>  R-squared:         </th> <td>   0.226</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>WLS</td>       <th>  Adj. R-squared:    </th> <td>   0.220</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   43.67</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Tue, 17 May 2022</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>11:14:05</td>     <th>  Log-Likelihood:    </th> <td> -19794.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  8906</td>      <th>  AIC:               </th> <td>3.971e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  8846</td>      <th>  BIC:               </th> <td>4.013e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>    59</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>    6.0258</td> <td>    0.152</td> <td>   39.556</td> <td> 0.000</td> <td>    5.727</td> <td>    6.324</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.36.0]</th>  <td>    0.0392</td> <td>    0.187</td> <td>    0.209</td> <td> 0.834</td> <td>   -0.328</td> <td>    0.406</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.37.0]</th>  <td>   -0.0815</td> <td>    0.189</td> <td>   -0.432</td> <td> 0.666</td> <td>   -0.452</td> <td>    0.289</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.38.0]</th>  <td>    0.0023</td> <td>    0.188</td> <td>    0.012</td> <td> 0.990</td> <td>   -0.366</td> <td>    0.371</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.39.0]</th>  <td>   -0.1169</td> <td>    0.188</td> <td>   -0.621</td> <td> 0.534</td> <td>   -0.486</td> <td>    0.252</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.40.0]</th>  <td>    0.2486</td> <td>    0.186</td> <td>    1.335</td> <td> 0.182</td> <td>   -0.116</td> <td>    0.614</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.41.0]</th>  <td>    0.1654</td> <td>    0.183</td> <td>    0.905</td> <td> 0.365</td> <td>   -0.193</td> <td>    0.524</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.42.0]</th>  <td>    0.2213</td> <td>    0.190</td> <td>    1.167</td> <td> 0.243</td> <td>   -0.151</td> <td>    0.593</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.43.0]</th>  <td>    0.1910</td> <td>    0.191</td> <td>    0.999</td> <td> 0.318</td> <td>   -0.184</td> <td>    0.566</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.44.0]</th>  <td>    0.3976</td> <td>    0.190</td> <td>    2.097</td> <td> 0.036</td> <td>    0.026</td> <td>    0.769</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.45.0]</th>  <td>    0.7522</td> <td>    0.187</td> <td>    4.027</td> <td> 0.000</td> <td>    0.386</td> <td>    1.118</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.46.0]</th>  <td>    0.6286</td> <td>    0.195</td> <td>    3.224</td> <td> 0.001</td> <td>    0.246</td> <td>    1.011</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.47.0]</th>  <td>    0.4826</td> <td>    0.190</td> <td>    2.540</td> <td> 0.011</td> <td>    0.110</td> <td>    0.855</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.48.0]</th>  <td>    0.5367</td> <td>    0.194</td> <td>    2.761</td> <td> 0.006</td> <td>    0.156</td> <td>    0.918</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.49.0]</th>  <td>    0.5856</td> <td>    0.195</td> <td>    3.008</td> <td> 0.003</td> <td>    0.204</td> <td>    0.967</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.50.0]</th>  <td>    0.5697</td> <td>    0.194</td> <td>    2.943</td> <td> 0.003</td> <td>    0.190</td> <td>    0.949</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.51.0]</th>  <td>    0.7385</td> <td>    0.198</td> <td>    3.727</td> <td> 0.000</td> <td>    0.350</td> <td>    1.127</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.52.0]</th>  <td>    0.6739</td> <td>    0.195</td> <td>    3.459</td> <td> 0.001</td> <td>    0.292</td> <td>    1.056</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.53.0]</th>  <td>    0.6313</td> <td>    0.202</td> <td>    3.125</td> <td> 0.002</td> <td>    0.235</td> <td>    1.027</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.54.0]</th>  <td>    0.7670</td> <td>    0.200</td> <td>    3.830</td> <td> 0.000</td> <td>    0.374</td> <td>    1.160</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.55.0]</th>  <td>    0.8264</td> <td>    0.204</td> <td>    4.055</td> <td> 0.000</td> <td>    0.427</td> <td>    1.226</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.56.0]</th>  <td>    0.7326</td> <td>    0.201</td> <td>    3.653</td> <td> 0.000</td> <td>    0.339</td> <td>    1.126</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.57.0]</th>  <td>    1.0893</td> <td>    0.204</td> <td>    5.329</td> <td> 0.000</td> <td>    0.689</td> <td>    1.490</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.58.0]</th>  <td>    0.9662</td> <td>    0.213</td> <td>    4.533</td> <td> 0.000</td> <td>    0.548</td> <td>    1.384</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.59.0]</th>  <td>    1.2559</td> <td>    0.224</td> <td>    5.596</td> <td> 0.000</td> <td>    0.816</td> <td>    1.696</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.60.0]</th>  <td>    0.8639</td> <td>    0.221</td> <td>    3.908</td> <td> 0.000</td> <td>    0.431</td> <td>    1.297</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.61.0]</th>  <td>    1.0555</td> <td>    0.223</td> <td>    4.741</td> <td> 0.000</td> <td>    0.619</td> <td>    1.492</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.62.0]</th>  <td>    1.3250</td> <td>    0.219</td> <td>    6.051</td> <td> 0.000</td> <td>    0.896</td> <td>    1.754</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.63.0]</th>  <td>    1.2049</td> <td>    0.229</td> <td>    5.263</td> <td> 0.000</td> <td>    0.756</td> <td>    1.654</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.64.0]</th>  <td>    1.2446</td> <td>    0.232</td> <td>    5.363</td> <td> 0.000</td> <td>    0.790</td> <td>    1.700</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.65.0]</th>  <td>    1.4242</td> <td>    0.234</td> <td>    6.084</td> <td> 0.000</td> <td>    0.965</td> <td>    1.883</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.66.0]</th>  <td>    0.9755</td> <td>    0.232</td> <td>    4.200</td> <td> 0.000</td> <td>    0.520</td> <td>    1.431</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.67.0]</th>  <td>    1.0702</td> <td>    0.243</td> <td>    4.399</td> <td> 0.000</td> <td>    0.593</td> <td>    1.547</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.68.0]</th>  <td>    1.2799</td> <td>    0.234</td> <td>    5.466</td> <td> 0.000</td> <td>    0.821</td> <td>    1.739</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.69.0]</th>  <td>    1.3089</td> <td>    0.245</td> <td>    5.345</td> <td> 0.000</td> <td>    0.829</td> <td>    1.789</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.70.0]</th>  <td>    1.5193</td> <td>    0.247</td> <td>    6.144</td> <td> 0.000</td> <td>    1.035</td> <td>    2.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.71.0]</th>  <td>    1.0835</td> <td>    0.253</td> <td>    4.281</td> <td> 0.000</td> <td>    0.587</td> <td>    1.580</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.72.0]</th>  <td>    1.3171</td> <td>    0.237</td> <td>    5.558</td> <td> 0.000</td> <td>    0.853</td> <td>    1.782</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.73.0]</th>  <td>    1.1630</td> <td>    0.237</td> <td>    4.901</td> <td> 0.000</td> <td>    0.698</td> <td>    1.628</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.74.0]</th>  <td>    1.4371</td> <td>    0.253</td> <td>    5.671</td> <td> 0.000</td> <td>    0.940</td> <td>    1.934</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(age)[T.75.0]</th>  <td>    1.3932</td> <td>    0.245</td> <td>    5.676</td> <td> 0.000</td> <td>    0.912</td> <td>    1.874</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(yr)[T.2004.0]</th> <td>   -0.0580</td> <td>    0.048</td> <td>   -1.218</td> <td> 0.223</td> <td>   -0.151</td> <td>    0.035</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_1</th>          <td>   -0.5309</td> <td>    0.127</td> <td>   -4.177</td> <td> 0.000</td> <td>   -0.780</td> <td>   -0.282</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_2</th>          <td>   -0.6819</td> <td>    0.096</td> <td>   -7.088</td> <td> 0.000</td> <td>   -0.870</td> <td>   -0.493</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_3</th>          <td>   -0.5519</td> <td>    0.087</td> <td>   -6.327</td> <td> 0.000</td> <td>   -0.723</td> <td>   -0.381</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_4</th>          <td>   -0.2330</td> <td>    0.084</td> <td>   -2.780</td> <td> 0.005</td> <td>   -0.397</td> <td>   -0.069</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_6</th>          <td>    0.2285</td> <td>    0.086</td> <td>    2.645</td> <td> 0.008</td> <td>    0.059</td> <td>    0.398</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_7</th>          <td>    0.2879</td> <td>    0.091</td> <td>    3.156</td> <td> 0.002</td> <td>    0.109</td> <td>    0.467</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_8</th>          <td>    0.4109</td> <td>    0.105</td> <td>    3.904</td> <td> 0.000</td> <td>    0.205</td> <td>    0.617</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_9</th>          <td>    0.4310</td> <td>    0.153</td> <td>    2.813</td> <td> 0.005</td> <td>    0.131</td> <td>    0.731</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ey10_10</th>         <td>    0.5187</td> <td>    0.147</td> <td>    3.540</td> <td> 0.000</td> <td>    0.231</td> <td>    0.806</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smoken</th>          <td>   -0.2906</td> <td>    0.068</td> <td>   -4.297</td> <td> 0.000</td> <td>   -0.423</td> <td>   -0.158</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokev</th>          <td>    0.2643</td> <td>    0.056</td> <td>    4.710</td> <td> 0.000</td> <td>    0.154</td> <td>    0.374</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>obese</th>           <td>    0.0901</td> <td>    0.054</td> <td>    1.662</td> <td> 0.097</td> <td>   -0.016</td> <td>    0.196</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>diabe</th>           <td>    0.7981</td> <td>    0.091</td> <td>    8.725</td> <td> 0.000</td> <td>    0.619</td> <td>    0.977</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>stroke</th>          <td>    0.8748</td> <td>    0.158</td> <td>    5.537</td> <td> 0.000</td> <td>    0.565</td> <td>    1.185</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hibpe</th>           <td>    1.0223</td> <td>    0.058</td> <td>   17.647</td> <td> 0.000</td> <td>    0.909</td> <td>    1.136</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>lunge</th>           <td>    0.7732</td> <td>    0.184</td> <td>    4.197</td> <td> 0.000</td> <td>    0.412</td> <td>    1.134</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hearte</th>          <td>    0.9840</td> <td>    0.083</td> <td>   11.867</td> <td> 0.000</td> <td>    0.821</td> <td>    1.147</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>male</th>            <td>   -0.9986</td> <td>    0.049</td> <td>  -20.546</td> <td> 0.000</td> <td>   -1.094</td> <td>   -0.903</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>1051.272</td> <th>  Durbin-Watson:     </th> <td>   1.805</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>  <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1484.205</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>           <td>-0.918</td>  <th>  Prob(JB):          </th> <td>    0.00</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>       <td> 3.792</td>  <th>  Cond. No.          </th> <td>    56.2</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            WLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                   logm   R-squared:                       0.226\n",
       "Model:                            WLS   Adj. R-squared:                  0.220\n",
       "Method:                 Least Squares   F-statistic:                     43.67\n",
       "Date:                Tue, 17 May 2022   Prob (F-statistic):               0.00\n",
       "Time:                        11:14:05   Log-Likelihood:                -19794.\n",
       "No. Observations:                8906   AIC:                         3.971e+04\n",
       "Df Residuals:                    8846   BIC:                         4.013e+04\n",
       "Df Model:                          59                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "===================================================================================\n",
       "                      coef    std err          t      P>|t|      [0.025      0.975]\n",
       "-----------------------------------------------------------------------------------\n",
       "Intercept           6.0258      0.152     39.556      0.000       5.727       6.324\n",
       "C(age)[T.36.0]      0.0392      0.187      0.209      0.834      -0.328       0.406\n",
       "C(age)[T.37.0]     -0.0815      0.189     -0.432      0.666      -0.452       0.289\n",
       "C(age)[T.38.0]      0.0023      0.188      0.012      0.990      -0.366       0.371\n",
       "C(age)[T.39.0]     -0.1169      0.188     -0.621      0.534      -0.486       0.252\n",
       "C(age)[T.40.0]      0.2486      0.186      1.335      0.182      -0.116       0.614\n",
       "C(age)[T.41.0]      0.1654      0.183      0.905      0.365      -0.193       0.524\n",
       "C(age)[T.42.0]      0.2213      0.190      1.167      0.243      -0.151       0.593\n",
       "C(age)[T.43.0]      0.1910      0.191      0.999      0.318      -0.184       0.566\n",
       "C(age)[T.44.0]      0.3976      0.190      2.097      0.036       0.026       0.769\n",
       "C(age)[T.45.0]      0.7522      0.187      4.027      0.000       0.386       1.118\n",
       "C(age)[T.46.0]      0.6286      0.195      3.224      0.001       0.246       1.011\n",
       "C(age)[T.47.0]      0.4826      0.190      2.540      0.011       0.110       0.855\n",
       "C(age)[T.48.0]      0.5367      0.194      2.761      0.006       0.156       0.918\n",
       "C(age)[T.49.0]      0.5856      0.195      3.008      0.003       0.204       0.967\n",
       "C(age)[T.50.0]      0.5697      0.194      2.943      0.003       0.190       0.949\n",
       "C(age)[T.51.0]      0.7385      0.198      3.727      0.000       0.350       1.127\n",
       "C(age)[T.52.0]      0.6739      0.195      3.459      0.001       0.292       1.056\n",
       "C(age)[T.53.0]      0.6313      0.202      3.125      0.002       0.235       1.027\n",
       "C(age)[T.54.0]      0.7670      0.200      3.830      0.000       0.374       1.160\n",
       "C(age)[T.55.0]      0.8264      0.204      4.055      0.000       0.427       1.226\n",
       "C(age)[T.56.0]      0.7326      0.201      3.653      0.000       0.339       1.126\n",
       "C(age)[T.57.0]      1.0893      0.204      5.329      0.000       0.689       1.490\n",
       "C(age)[T.58.0]      0.9662      0.213      4.533      0.000       0.548       1.384\n",
       "C(age)[T.59.0]      1.2559      0.224      5.596      0.000       0.816       1.696\n",
       "C(age)[T.60.0]      0.8639      0.221      3.908      0.000       0.431       1.297\n",
       "C(age)[T.61.0]      1.0555      0.223      4.741      0.000       0.619       1.492\n",
       "C(age)[T.62.0]      1.3250      0.219      6.051      0.000       0.896       1.754\n",
       "C(age)[T.63.0]      1.2049      0.229      5.263      0.000       0.756       1.654\n",
       "C(age)[T.64.0]      1.2446      0.232      5.363      0.000       0.790       1.700\n",
       "C(age)[T.65.0]      1.4242      0.234      6.084      0.000       0.965       1.883\n",
       "C(age)[T.66.0]      0.9755      0.232      4.200      0.000       0.520       1.431\n",
       "C(age)[T.67.0]      1.0702      0.243      4.399      0.000       0.593       1.547\n",
       "C(age)[T.68.0]      1.2799      0.234      5.466      0.000       0.821       1.739\n",
       "C(age)[T.69.0]      1.3089      0.245      5.345      0.000       0.829       1.789\n",
       "C(age)[T.70.0]      1.5193      0.247      6.144      0.000       1.035       2.004\n",
       "C(age)[T.71.0]      1.0835      0.253      4.281      0.000       0.587       1.580\n",
       "C(age)[T.72.0]      1.3171      0.237      5.558      0.000       0.853       1.782\n",
       "C(age)[T.73.0]      1.1630      0.237      4.901      0.000       0.698       1.628\n",
       "C(age)[T.74.0]      1.4371      0.253      5.671      0.000       0.940       1.934\n",
       "C(age)[T.75.0]      1.3932      0.245      5.676      0.000       0.912       1.874\n",
       "C(yr)[T.2004.0]    -0.0580      0.048     -1.218      0.223      -0.151       0.035\n",
       "ey10_1             -0.5309      0.127     -4.177      0.000      -0.780      -0.282\n",
       "ey10_2             -0.6819      0.096     -7.088      0.000      -0.870      -0.493\n",
       "ey10_3             -0.5519      0.087     -6.327      0.000      -0.723      -0.381\n",
       "ey10_4             -0.2330      0.084     -2.780      0.005      -0.397      -0.069\n",
       "ey10_6              0.2285      0.086      2.645      0.008       0.059       0.398\n",
       "ey10_7              0.2879      0.091      3.156      0.002       0.109       0.467\n",
       "ey10_8              0.4109      0.105      3.904      0.000       0.205       0.617\n",
       "ey10_9              0.4310      0.153      2.813      0.005       0.131       0.731\n",
       "ey10_10             0.5187      0.147      3.540      0.000       0.231       0.806\n",
       "smoken             -0.2906      0.068     -4.297      0.000      -0.423      -0.158\n",
       "smokev              0.2643      0.056      4.710      0.000       0.154       0.374\n",
       "obese               0.0901      0.054      1.662      0.097      -0.016       0.196\n",
       "diabe               0.7981      0.091      8.725      0.000       0.619       0.977\n",
       "stroke              0.8748      0.158      5.537      0.000       0.565       1.185\n",
       "hibpe               1.0223      0.058     17.647      0.000       0.909       1.136\n",
       "lunge               0.7732      0.184      4.197      0.000       0.412       1.134\n",
       "hearte              0.9840      0.083     11.867      0.000       0.821       1.147\n",
       "male               -0.9986      0.049    -20.546      0.000      -1.094      -0.903\n",
       "==============================================================================\n",
       "Omnibus:                     1051.272   Durbin-Watson:                   1.805\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1484.205\n",
       "Skew:                          -0.918   Prob(JB):                         0.00\n",
       "Kurtosis:                       3.792   Cond. No.                         56.2\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = wls('logm ~ ey10_1 +ey10_2 + ey10_3 + ey10_4 + ey10_6 + ey10_7 + ey10_8 + ey10_9 + ey10_10 + C(age) + C(yr)  + smoken + smokev + obese + diabe + stroke + hibpe + lunge + hearte + male',data=data,var_weights=np.asarray(data['perwt'])).fit()\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    -0.530896\n",
       "2    -0.681907\n",
       "3    -0.551857\n",
       "4    -0.233032\n",
       "5     0.000000\n",
       "6     0.228472\n",
       "7     0.287888\n",
       "8     0.410882\n",
       "9     0.431045\n",
       "10    0.518710\n",
       "dtype: float64"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_m_data_h = model.params[['ey10_'+str(q) for q in range(1,11) if q!=5]]\n",
    "grad_m_data_h.index = [x for x in range(1,11) if x!=5]\n",
    "grad_m_data_h[5] = 0\n",
    "grad_m_data_h = grad_m_data_h.sort_index()\n",
    "grad_m_data_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = pd.DataFrame(index=[x for x in range(1,11)],columns=['sim','mean','se'])\n",
    "table['sim'] = grad_m_sim \n",
    "table['mean_h'] = grad_m_data_h\n",
    "table['mean'] = grad_m_data\n",
    "table['se'] = grad_m_data_se "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sim</th>\n",
       "      <th>mean</th>\n",
       "      <th>se</th>\n",
       "      <th>mean_h</th>\n",
       "      <th>up</th>\n",
       "      <th>low</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-0.739822</td>\n",
       "      <td>-0.310938</td>\n",
       "      <td>0.156416</td>\n",
       "      <td>-0.530896</td>\n",
       "      <td>-0.004362</td>\n",
       "      <td>-0.617514</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.663559</td>\n",
       "      <td>-0.524490</td>\n",
       "      <td>0.114373</td>\n",
       "      <td>-0.681907</td>\n",
       "      <td>-0.300318</td>\n",
       "      <td>-0.748661</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.559215</td>\n",
       "      <td>-0.429465</td>\n",
       "      <td>0.095196</td>\n",
       "      <td>-0.551857</td>\n",
       "      <td>-0.242880</td>\n",
       "      <td>-0.616049</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-0.366673</td>\n",
       "      <td>-0.218555</td>\n",
       "      <td>0.090118</td>\n",
       "      <td>-0.233032</td>\n",
       "      <td>-0.041925</td>\n",
       "      <td>-0.395186</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.055878</td>\n",
       "      <td>0.182247</td>\n",
       "      <td>0.085721</td>\n",
       "      <td>0.228472</td>\n",
       "      <td>0.350260</td>\n",
       "      <td>0.014234</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>-0.032321</td>\n",
       "      <td>0.193383</td>\n",
       "      <td>0.090745</td>\n",
       "      <td>0.287888</td>\n",
       "      <td>0.371243</td>\n",
       "      <td>0.015524</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>-0.017302</td>\n",
       "      <td>0.213062</td>\n",
       "      <td>0.100431</td>\n",
       "      <td>0.410882</td>\n",
       "      <td>0.409907</td>\n",
       "      <td>0.016217</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.036369</td>\n",
       "      <td>0.136299</td>\n",
       "      <td>0.127536</td>\n",
       "      <td>0.431045</td>\n",
       "      <td>0.386269</td>\n",
       "      <td>-0.113671</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.121999</td>\n",
       "      <td>0.187443</td>\n",
       "      <td>0.131535</td>\n",
       "      <td>0.518710</td>\n",
       "      <td>0.445253</td>\n",
       "      <td>-0.070366</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         sim      mean        se    mean_h        up       low\n",
       "1  -0.739822 -0.310938  0.156416 -0.530896 -0.004362 -0.617514\n",
       "2  -0.663559 -0.524490  0.114373 -0.681907 -0.300318 -0.748661\n",
       "3  -0.559215 -0.429465  0.095196 -0.551857 -0.242880 -0.616049\n",
       "4  -0.366673 -0.218555  0.090118 -0.233032 -0.041925 -0.395186\n",
       "5   0.000000  0.000000  0.000000  0.000000  0.000000  0.000000\n",
       "6   0.055878  0.182247  0.085721  0.228472  0.350260  0.014234\n",
       "7  -0.032321  0.193383  0.090745  0.287888  0.371243  0.015524\n",
       "8  -0.017302  0.213062  0.100431  0.410882  0.409907  0.016217\n",
       "9   0.036369  0.136299  0.127536  0.431045  0.386269 -0.113671\n",
       "10  0.121999  0.187443  0.131535  0.518710  0.445253 -0.070366"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table['up'] = table['mean'] + 1.96*table['se']\n",
    "table['low'] = table['mean'] - 1.96*table['se']\n",
    "table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/h9/sn9_m3953j93s5pv8q78dfvh0000gn/T/ipykernel_37762/3838045957.py:5: UserWarning: linestyle is redundantly defined by the 'linestyle' keyword argument and the fmt string \"-D\" (-> linestyle='-'). The keyword argument will take precedence.\n",
      "  ax.plot(table.index,table['mean_h'],'-D',color='r',linestyle=':',label='Data (controls for health)')\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABq3ElEQVR4nO2dd3iUxfbHP5NNbxAICSV0gat0pCoKCHJVEFH0KiBdsaGigu13FcvVq4JYrl4URUQEQb0IiL2goAYUFFFRQCFASCWkZ0uye35/bCFlN9kku2nM53n22bfMO3N2s3m/78ycOUeJCBqNRqPReCKgvg3QaDQaTcNGC4VGo9FoKkULhUaj0WgqRQuFRqPRaCpFC4VGo9FoKiWwvg3wB7GxsdKpU6f6NkOj0WgaDbt37z4hIq3cnWuSQtGpUyd27dpV32ZoNBpNo0EpdcTTOT30pNFoNJpK0UKh0Wg0mkrRQqHRaDSaSmmScxTuKC4uJjk5GZPJVN+maE5TQkNDSUhIICgoqL5N0WiqxWkjFMnJyURFRdGpUyeUUvVtjuY0Q0TIysoiOTmZzp0717c5mqbI1q0waxasXAmjRvm06tNm6MlkMtGyZUstEpp6QSlFy5YtdY9W4x+2boXx4+HIEfv71q0+rf60EQpAi4SmXtG/P41fcIpEUZF9v6jI52JxWgmFRqPRNCnKi4QTH4uFFgqNRqNprMyaVVEknBQV2c/7AC0UHlizBjp1goAA+/uaNbWv02Aw0K9fP3r27Enfvn1ZunQpNput0muSkpJYu3ZttdsyGo2MGDECq9VaU3N9zrPPPkuRpx91JURGRlZ6fsyYMWRnZ9fULI2m8fHJJ7Bxo33iOjzcfZnwcPt5H6CFwg1r1sDcufZ5IRH7+9y5tReLsLAw9uzZw2+//cZnn33Ghx9+yMMPP1zpNTUVitdee40rrrgCg8FQU3N9TmVCURtBmzZtGv/9739rfL1G0ygwGk9t//vf8OSTdu+mLVsqikV4uP24r7yfRKTJvc4++2wpz759+1zbt98uMmKE51dIiIhdIsq+QkI8X3P77RWarEBERESZ/b/++ktatGghNptNDh8+LMOHD5f+/ftL//795dtvvxURkSFDhkh0dLT07dtXli5d6rFceYYNGyaHDx8WEZGtW7fKiBEjZNKkSdKjRw+ZMmWK2Gw2ERH5/PPPpV+/ftKrVy+ZNWuWmEymCnUdPHhQRo8eLX369JH+/fvLn3/+KTabTRYsWCA9e/aUXr16ybp16ypt67nnnpOgoCDp1auXjBw50vV9PPDAAzJ48GDZvn27PP3009KzZ0/p2bOnPPPMMxW+t5SUFDnvvPOkb9++0rNnT9m2bZuIiJw8eVJ69uxZ9R+gAVD6d6jReM0LL4i0bCmSn2/fP3pUxGw+df7LL0XCw+03qvBw+341AXaJh3tqvd/U/fGqrVC4Ewnny5dCISLSvHlzSUtLk8LCQjEajSIicuDAAXF+hq1bt8q4ceNc5T2VK43ZbJb4+HjX/tatWyU6OlqOHTsmVqtVhg4dKtu3bxej0SgJCQmyf/9+ERGZNm1amRu0k8GDB8uGDRtERMRoNEphYaG8++67MmbMGCkpKZG0tDRp3769pKSkeGxLRKRjx46SmZnpqheQ9evXi4jIrl27pFevXlJQUCD5+fly1llnyY8//ljme1uyZIn861//EhGRkpISycvLc9V1xhlnyIkTJzx/+Q0ELRQar8jOFlmyROTIEfv+jh0id94pUtlv/MsvRTp2rJFIiFQuFKfNgrvSPPts5ec7dbIPN5WnY0f46ivf2mL/+9hXjs+bN489e/ZgMBg4cOCA2/LelDtx4gTNmzcvc2zw4MEkJCQA0K9fP5KSkoiKiqJz5850794dgBkzZvDiiy8yf/5813X5+fkcP36cyy+/HLCvLgb45ptvmDx5MgaDgfj4eEaMGMEPP/xAdHS027aGDx9ewU6DwcCkSZNc9V1++eVEREQAcMUVV7B9+3b69+/vKj9o0CBmz55NcXExEydOpF+/fq5zcXFxpKSk0LJlS/dftEbT0BGBwkKIjITcXFi40L59ww0wZIj9VRmjRkFSkl9M03MUbnjsMfdDfo895tt2Dh06hMFgIC4ujmeeeYb4+Hh+/vlndu3ahcVicXuNN+XCwsIqLOwKCQlxbRsMBkpKSlwiVRmeylR2rbu23BEaGuqaQ/HGlvPPP59t27bRrl07pk2bxhtvvOE6ZzKZCAsLq7IOjaZBIgLDhsG8efb9jh3h8GG7SDQAtFC4YepUWL7c/rdSyv6+fLn9uK/IzMzkxhtvZN68eSilyM3NpU2bNgQEBLB69WrX5G5UVBT5+fmu6zyVK01MTAxWq7XKVcB/+9vfSEpK4s8//wRg9erVjBgxokyZ6OhoEhIS2LhxIwBms5mioiLOP/981q9fj9VqJTMzk23btjF48OBK2yv/WUpz/vnns3HjRoqKiigsLOS9997jvPPOK1PmyJEjxMXFcf311zNnzhx+/PFHwC4yaWlp6GRVmkbFjz/aJ6TBfqO56iq48MJT5zt2rB+73HBaDj15w9SpvhUGsLus9uvXj+LiYgIDA5k2bRp33nknADfffDOTJk3inXfeYdSoUa4hmD59+hAYGEjfvn2ZOXOmx3LlGTt2LN988w1jxozxaE9oaCgrV67kqquuoqSkhEGDBnHjjTdWKLd69WpuuOEGHnzwQYKCgnjnnXe4/PLLSUxMpG/fviileOqpp2jdujV//PGHx/bmzp3LxRdfTJs2bdhabiHQgAEDmDlzpktsrrvuujLDTgBfffUVixcvJigoiMjISFePYvfu3QwdOpTAQP1z1jRw8vMhIsLud//JJ/D44zBnDsTGwl131bd1HlHedPkbGwMHDpTyGe5+//13zjzzzHqyqO756aefWLp0KatXr65vU/zO7bffzoQJExg9enR9m1Ilp9vvUFOKnTth9GjYtMn+npdnPx4dXb92OVBK7RaRge7O6aGnJkr//v0ZNWpUg1pw5y969erVKERC00TZutXuAVM+XIbVCitW2NczAPTtC9OnQ9u29v3o6AYjElWhexQaTR2if4dNjNKxlpyL3Pr1g5gY+wR1z54wYAC8+abfTbHZbBQXF5dxJqkOlfUo9KCuRqPR1AR3UVvHjoXmzSE1FQID7f70rVr5pXkRobi4mKKiIvLy8lyOIt27d/d5cqxqDT0ppWKUUn18aoFGo9HUFE/DPr7EYgGni3dKin04acMG91FbS0rsayC++MK+Hxdn92jyESUlJRQUFJCWlsb+/fvZv38/x44do6ioiLCwMK/czGtClUKhlPpKKRWtlGoB/AysVEot9Ys1Go1G4y2+SNZTWAiffgrHj9v3Dx+2uzs6h66//RZCQk6ttD1wAK67Dm6+2XPU1uJin61/sNlsGI1GTpw4wV9//cUff/zB4cOHOXnyJEFBQURFRREVFUVoaCgBAf6bcvam5mYikgdcAawUkbMBzz6XTYitW7fSqVOnCq6cNUUpxbRp01z7JSUltGrVivHjx1ernk6dOnHixIlal9FoGi2ekvV88QX88Yd96AegoADuuw+2bbPvp6TYJ5NXrbLvp6fD3/8On31m3xeBHTsgI8O+37UrPPzwqTUNQ4bYhWn1ar9EbRURzGYzOTk5HD16lN9//50///yTtLQ0bDYbERERREVFER4eXqcBP70RikClVBvgH8AWXzaulLpIKbVfKfWnUureSsoNUkpZlVJX+rL9yti6dSvjx4/nyJEjjB8/3idiERERwa+//orREQXys88+o127drWuV6M5ragsWc+ECXDmmfDqq/ZjQUHw9NP2xW0ALVrAuHH24SqAhATYvh0uvdS+36UL/PUXXHKJfb91a3jwQejWzb4fFgYdOtgXxvkoaqtzOCk1NZX9+/dz4MABkpOTXcNJUVFRREZGEhQUVG9ZEr0RikeAT4C/ROQHpVQX4GBtG1ZKGYAXgYuBs4DJSqmzPJR70mFDneAUCWdI7KKiIp+JxcUXX8wHH3wAwFtvvcXkyZNd506ePMnEiRPp06cPQ4cOZe/evQBkZWUxduxY+vfvzw033FBmHPLNN99k8ODB9OvXjxtuuOG0cIfVnL6ICHLttZUm65GWLbFcdhlWqxUJDgazGZzxy0JD4ZVXwBmBIDgYhg+HmsQIKx/i20uRsNlsFBUVVRhOys7Odg0nRUZG+n04qTpU6fUkIu8A75TaPwRM8kHbg4E/HfWhlFoHXAbsK1fuVuB/wCAftAnA/Pnz2bNnj9tz2dnZ/PrrrxUSChUVFTFmzBh69epFTExMhev69evHs1VFGwSuueYaHnnkEcaPH8/evXuZPXs227dvB2DRokX079+fjRs38uWXXzJ9+nT27NnDww8/zPDhw3nwwQf54IMPWL58OWB3tVy/fj3ffvstQUFB3HzzzaxZs4bp06dX7wvRaBogVqsVi8WCxWLBaDRSVFSE0WikQ6dORKak4O7Z2hYaypGnnqIwIAB+/x2lFEFBQYSEhJR5NxgMBAYGYjAYMBgMNX9Sd4rFrFn24SY3IiEiWCwWl3dSQUEBNpsNpRQhISFEREQ0+HzqVQqFUqo7sAyIF5FeDq+nCSLyr1q23Q44Vmo/GSgTHlEp1Q64HLiAKoRCKTUXmAvQoUOHGhu1f/9+j1nnbDYb+/fvZ+jQoTWuv0+fPiQlJfHWW29xibN76+Cbb77hf//7HwAXXHABWVlZ5Obmsm3bNjZs2ADAuHHjXEL1xRdfsHv3bgYNsn81RqORuLi4Gtum0dQHTv9/i8WC2WymsLAQo9F4Kpik1UqLTZtQ551HWMeO5Dz7LMadO4m9+24CSiXzsYWFkbFiBWrYMJw5EUUEq9XqilFmtVoRkQo35uDgYIKCgggODnaJiVNInO8eb+ZuoraWlJRgNBrJz88nLy/P9VkCAwMJCwtrMD0Fb/FmHcUrwELgZQAR2auUWgvUVijcfevlfbueBe4REWtViisiy4HlYF9wV1nZyp78yw87lSY8PJwtW7YwqpZZoyZMmMCCBQv46quvyMrKch1359rm/NzuPr+IMGPGDP7973/Xyh6Npi4QEUpKSiguLnbduI1GIyaTyXXzVkoRGBhIUFCQK6S9IS2N1v/+N3lZWeTccQfWVq0oGj+ejJYtiZszhwCj0SUSpmHDyrTprK8qu5xiZTKZsNlsrid+p10i4hKTkJAQgoODCQ4OLiMkxcXFFBYWkpub6/pMgYGBBAcHuz5LY8UboQgXke/L3ajcx42uHslA+1L7CUBKuTIDgXWOtmOBS5RSJSKy0Qftu2XUqFFs2bKlglj4SiQAZs+eTbNmzejduzdflUpwcf7557NmzRoeeOABvvrqK2JjY4mOjnYd/+c//8lHH33kyg89evRoLrvsMu644w7i4uI4efIk+fn5dGxAUSc1pydWq9XVSzCZTK5eQuneuvMm627oJfDwYcK3biVv9mysrVuTumkTxc4JZQemYcPIWLGC2IULObF4cQWR8BallGsIyhNOMSkpKcFisWC1Wl2fpbztjWU4qTp4IxQnlFJdcTztOzyPUn3Q9g9AN6VUZ+A4cA0wpXQBEens3FZKvQ5s8adIOCkvFr4UCYCEhARuv/32CscfeughZs2aRZ8+fQgPD2eVw4Vv0aJFTJ48mQEDBjBixAjX0NpZZ53Fv/71L8aOHYvNZiMoKIgXX3xRC4WmznCuDnbOJTgFwWw2l+kNO3sI3g65RG7cSPRrr1Fw6aXYWrWi2JFcqzymYcNI/uYbn30eT3gjJk2ZKmM9ObyclgPnANnAYWCqiLjJAVfNxpW6BPvwkgF4TUQeU0rdCCAiL5Ur+zp2oXi3qnp9Fetp69atzJo1i5UrV/pMJDSnN40p1pNzfL/0yzl05BymcQ7VOEXBOWxU7RtqSQlRa9di7t0bS//+qKIiAgoLsfop/EVTpaCggB49etQohEeNYz05XFNvEpExSqkIIEBE3GeeqQEi8iHwYbljL3koO9NX7XrLqFGjSPJTakGNpr6w2WxuBcDZK3D2EEpKSioMn4gIBoOBgIAADAYD4eHhPhliUWYzzV58kaJx4zjZvz8SHo7V04I2TZ1TqVA4JpHPdmwX1o1JGo2mujjH0MsLgPOmX/rdarWWubk7J2ydAhAQEOCagPXnOHvQwYNEvv022fffj0REkLppE9b4eL+1p6k53sxR/KSU2ox9LYVLLERkg9+s0mg0bhERcnNzKzz9FxcXuy3rvPE7RaAhLeIK2bOHqPXryZ88mZIuXbC2bl3fJmk84I1QtACysK9lcCKAFgqNpo45efIkycnJLpfMgIAAlwtmg/eyKSkheuVKStq3p+iiiyiYNImiMWOwuVnAqmlYeLMye1ZdGKLRaCqnuLiYtLQ0IiMjG6f3jVJEbNqEpVcvii66CAICtEg0ErxZmb2SigvhEJHZfrFIo9G4JS0tDaBRiUTQ/v00W7aMrCeeQEJDSVu7Fmkk6T81p/BmsHIL8IHj9QUQDRT406gGQ10kRdFovKCgoIDs7GzCG5knkCE7m7Bt2wg6cABAi0QjpUqhEJH/lXqtwR5uvJf/TatnfJEUpRwGg4F+/frRs2dP+vbty9KlSz3GlXKSlJTE2rVrq92W0WhkxIgRPo0mu2fPHj788MOqC5bjoYceYsmSJV6Xnzx5Mn369OGZZ56pdls1aa8yRo4ciXNNzuOPP+46npSURK9e7v8NFixYwJdffumT9sG+yvn48eN+90LyCSI0e/FFol97DQDT0KEkf/MNlj46MWZjpibuD92Amkfdawx4SopSS7EICwtjz549/Pbbb3z22Wd8+OGHPPzww5VeU1OheO2117jiiit8OkxRmVC4ArjVkrS0NL777jv27t3LHXfc4dU1vmq7KkoLRWXceuutPPHEEz5rNysrC4vFQnBwsM/qrC2hiYkkDB9OaGJi2RNKEbx3L8G//eY6JI2sF6SpiDepUPOVUnnOF/A+cI//TfMzI0fC66/bt4uL7ftvvmkXg3Hj3CdFufhi+/kTJ+zl33/ffs4xdlwd4uLiWL58OS+88AIiQlJSEueddx4DBgxgwIABfPfddwDce++9bN++nX79+vHMM894LFeeNWvWcNlll7n2n3rqKXr37k3fvn259157jqg9e/YwdOhQ+vTpw+WXX+6KITVy5EjuueceBg8eTPfu3dm+fTsWi4UHH3yQ9evX069fP9avX89DDz3E3LlzGTt2LNOnT+fIkSOMHj2aPn36MHr0aI4ePVrBrueff56zzjqLPn36cM0111Q4P3bsWDIyMujXrx/bt2+v1Mb777+fESNG8Nxzz1WoZ9++fYwcOZIuXbrw/PPPu457yt9x0003MXDgQHr27MmiRYsq1HfvvfdiNBrp168fU6dOBexP+tdffz09e/Zk7NixroRUHTt2JCsryzWnUBtMJhPp6elERETUui5fEZqYSNycOQQeP07cnDlEvP028dOnY8jMBCDz+ec58fTT9Wylxpd4M/QUJSLRpV7dReR/dWFcvTBrFpQKXVwGs9l+3kd06dIFm81GRkYGcXFxfPbZZ/z444+sX7+e2267DYAnnniC8847jz179riC/7krVxqLxcKhQ4fo5Mji9dFHH7Fx40Z27tzJzz//zN133w3A9OnTefLJJ9m7dy+9e/cu07spKSnh+++/59lnn+Xhhx8mODiYRx55hKuvvpo9e/Zw9dVXA7B79242bdrE2rVrmTdvHtOnT2fv3r1MnTrVrW1PPPEEP/30E3v37uWllyouwt+8eTNdu3Zlz549nHfeeZXamJOTw9dff81dd91VoZ4//viDTz75hO+//56HH36Y4uLiMvk79uzZg8FgYM2aNQA89thj7Nq1i7179/L111+7kkaVttvZI3Rec/DgQW655RZ+++03mjdv7goRDzBgwAC+/fZbN3917xERUlNTCQwMbDBrH5wi4QzvHWA00vLBBwnet4/Aw4fthUJC6tFCjT/wxuvpCxEZXdWxRkepqK0EBZ3ab9fOfZpFOJULNza27PW1WCjkjLVVXFzMvHnzXDewA47Jv/J4U+7EiRM0b97ctf/5558za9Ys10RoixYtyM3NJScnhxGOTF8zZszgqquucl1zxRVXAHD22WdXGsZkwoQJhIWFAZCYmOjKmzFt2jSXIJWmT58+TJ06lYkTJzJx4kSP9QJV2ugUK3eMGzeOkJAQQkJCiIuLIz09vdL8HW+//TbLly+npKSE1NRU9u3bR58qxtU7d+5Mv379gIrfU1xcHCkp5YMhV4/c3Fzy8/OJbiATwOVFwkmA2YxNKZTOrthk8SgUSqlQIByIVUrFcCp/RDTQtg5sqx+cGavKi0UNcuFWxaFDhzAYDMTFxfHwww8THx/Pzz//jM1m8xi//plnnqmyXFhYGCaTybXvLlFLVYQ4ngoNBkOlcwCVDYm4a/ODDz5g27ZtbN68mUcffZTffvutynwBNWk7pNRTrfMzeMrfcfjwYZYsWcIPP/xATEwMM2fOLPP9eduGsdQN1GQyuQS0JpSUlJCSktJgvJxCt28n7qabKoiEkwCTidiFC+skkqum7qmsP3sDsBv4G/CjY3s3sAl7ruumSw1z4VaHzMxMbrzxRubNm4dSitzcXNq0aUNAQACrV692jZ1HRUWRn38qDqOncqWJiYnBarW6bnZjx47ltddec+XXOHnyJM2aNSMmJsaVhnX16tWuJ3dPlLelPOeccw7r1q0D7HMkw4cPL3PeZrNx7NgxRo0axVNPPUVOTg4FBZ49rWtiY2WMHj2ad999l4yMDMD+PRw5coS8vDwiIiJo1qwZ6enpfPTRR26vDwoKchsqwx0HDhzw6BXlDenp6a7EN/VB0IEDxDz6KMrxGwo+eBAJCsLm4QHGFhbGicWL69JETR3iUShE5DlHPogFItK51KuviLxQhzbWD06x6NjRZyLhnAzt2bMnY8aMYezYsa6J05tvvplVq1YxdOhQDhw44Hpa7tOnD4GBgfTt25dnnnnGY7nyjB07lm8cT3cXXXQREyZMYODAgfTr18/lOrpq1SoWLlxInz592LNnDw8++GAVX8ko9u3b55rMLs/zzz/PypUr6dOnD6tXr64wyWy1Wrn22mvp3bs3/fv354477igzROaO6tpYGaXzd/Tp04cLL7yQ1NRU+vbtS//+/enZsyezZ8/m3HPPdXv93LlzXUNnlVFcXMyff/7JwIFuIzZXSVFREVlZWXXamwjIzKTZsmUEJicDEJiSQtSaNQQdPAhA3rXXcuzHH8l47TVs5XpKnrLLaZoOHvNRKKUuEJEvlVJXuDvfkIMC+iofRWPmp59+YunSpaxevbq+TTnteO+99/jxxx959NFHK5yr6ndos9n466+/XKk3/YUymwnfsoXi7t2x9O5N4NGjJIwYQebTT1N4xRVQXIyy2RA3E9Ol5yq0SDQs/JWPorKhJ2cf/1I3r/HVtkJTp/Tv359Ro0b5dMGdxjtKSkrcemJ5w8mTJzGZTH4RidDt2wndsQOwx+Rp+cADRDhcvEs6dODo99/bRQIgKMitSMCpFKQl7dppkThNqDLDXWPEU4/ib3/7W8Nf2appsogIf/zxh8cehcVi4cCBA4SHh/vEHTbo4EEMKSmYHPM6bf/+d6ytW5PuSLEbePQoJQkJ0EBcbzW1p84z3Cml7qysUhFZWm1L6pHQ0FCysrJo2bKlFgtNnSMiZGVlefRmc66ZcOaP8ERoYiKxCxdyYvHiCk/yASdOEPLbbxgdwtB86VJC9u61eyIpRcZ//4u1TRtX+ZIOTTvAgsZ3VOZSEeV47wEMAjY79i8FtvnTKH+QkJBAcnIymY7VoxpNXRMaGkpCQoLbc/n5+eTm5hIVFeX2PJSdG4ibM4eMZcsgOBjTkCEQEED0ypU0W76coz/+iERFkb1ggT18huPBqKRrV798Lk3Tp8qhJ6XUp8AkZ65spVQU8I6IXFQH9tUId0NPGk1DxWq1cuDAAQIDAz0OGbhb7CZBQajiYlI++ADLWWcReOwYAdnZWHr10sNJpyn1MZntpANgKbVvATpV2wqNRuOWjIwMrFZrtUQCQBUXI8HBBDjWhZS0b2+P0nqaiMTGjREMH55Aly4dGT48gY0bG048rKaGN6t5VgPfK6Xew+4scTnwhl+t0mhOE4xGI5mZmURGRnosE7twoccV0cpiIfaf/zztVkRv3BjB/fe3xGi0i+Lx44Hcf39LACZOLKxP05ok3qRCfUwp9THgXGY7S0R+8kXjSqmLgOcAA/CqiDxR7vxUTkWqLQBuEpGffdG2RlPfiAjHjx8nODi40gns7IULaXnvvQS4CSvSlFdEm82QlWXgxAkDmZn2d+dr/fpIl0g4MRoDWLIkRguFH/A2PsAeINVZXinVQUQqxpCuBkopA/ZQIBcCycAPSqnNIrKvVLHDwAgRyVZKXQwsB4bUpl2NpqGQnZ2N0WisdAKbkhKaP/ccxV27EnToUJmeRX0tdtu4MYIlS2JISTHQtq2VBQuyvb45V3bzP7UfQGamgbw897lUoqJsGI3uPRePHzeQnBxIQkLd5Cg5XfAmeuytwCIgHbBiDw4oQG1TVg0G/hSRQ4521gGXAS6hEJHSyRZ2AO5dRjSaRkZxcTGpqalVh+kIDCTr4YexxcQQkJ9f7yui3Q353HdfS3JyAujf3+zhxu/dzT821kpsrJXu3Ys55xwTrVpZXcdiY62OfRuhocLw4QkcP+7u9qU477wEBg40MWFCIePGFdKiReVZJJsCpcW7fXt4/HGoItJMtfDG6+lPYIiIZPmuWVBKXQlcJCLXOfanOdqZ56H8AuBvzvJuzs8F5gJ06NDh7CNHjvjSXI3Gpxw7dswVjNAdATk5BP3xB+ahQ13HNm6M4LvHfuapE9dxd+yrnPN/fWs1zCICJpOioEBRWBhAQUGAa9u+X3a7oCCAzZsjKgz5eKL0zf/Ujd7dvv3mXx3KCxZAWJiNBQuyMZkC2LQpggMHggkMFM4/38iECYVceGER4eFNb4Gxu+8iPByWL6+eWFTm9eSNUGwFLhQRn/bllFJXAX8vJxSDReRWN2VHAf8FhnsjWNo9VtOQKSgo4NChQ0RFRXlc/NnynnuI+PBDkrdvx9a8udubQUiIjZtuymXAALPjhq4cN/tT26WPl775FxQEUFSksFq9W3waHm4jIsJGZqaBUxkHSiO88kpGKRGwERLi35tyZUNgIvDHH0Fs2hTJ++9HkJISSHi4jQsvLOKyywoZPtxIDTxIGxxms2L48HacOFGxd9WxI1SSSqYCtRWKFdgX3X0AmJ3Ha7syWyk1DHhIRP7u2L/PUe+/y5XrA7wHXCwi7rP5lEMLhaahYrVa+fPPP1FKVRrPKSA3l+B9+1xDS56HWtwTGChERtqIiHC+27cjImxERtqIjDy1fep4xfMRETbCwwVn6nVPdrRrV8I33yRX78uoI2w2+OGHEDZtiuTDD8PJzTXQooWVceMKueyyQgYMMNNYgjXk5gawe3cIu3aF8MMPoezdG4LF4t54peyf3VtqKxQVEwgDIvKwu+PVMCoQOACMBo4DPwBTROS3UmU6AF8C08vNV1SKFgpNQyUjI4P09HT3E9gWC1Fr1pA/fTquO7ODzp074ulJ/u2308rd6IXgYPHLzc/TkM/jj2c1Cm8jiwW+/jqMTZsi+fzzMMzmABISipkwoZCJEwvp1s27fCN1RUqKgV27QvnhB7swHDgQhIgiMFDo3dvMwIFm/ve/SE6erDj348sehTfusQ87KokQEZ/9EkSkRCk1D/gEu3vsayLym1LqRsf5l4AHgZbAfx1d9BJPH0SjaeiYTCYyMjI8zkuEf/opLR95hOJu3TCVSvr0xhuevaLatbMyaJDZ43lf4xSDmno91TfBwXDhhUYuvNBIQYHik0/C2bw5kpdeasZ//9ucs84yM2FCIRMmFNKmTd1GXrbZ4M8/g/jhh1PCkJJiv0VHRNgYMMDMuHGFDBpkpm9fM2Fh9of8s86yuJ2jeOwx39nmTY9iGLACiBSRDkqpvsANInKz78zwLbpHoWloiAhJSUlVpkgN/uUXLL17A1BSAo880oLVq6M56ywzhw4FYTI1zif5hk5mZgBbtkSweXMke/aEoJQwZIjdc+qSS4po1sz3nlMWC/zySwg//BDCrl2h7N4dQk6OvWfQqlUJgwaZGTjQxODBZnr0sFBZssOKXk+q2l5PtR162glcCWwWkf6OY7+KSM3zPPoZLRSahkZOTg5Hjx4lOjq67AkRmj3/PIUTJ1LSsaPrcF6eYt68OLZvD+P663O5555s3n+/5usXNN6TlBTI5s0RbNwYyeHDQQQHCyNHFjFhQiGjRxur7aHlJC9P8dNPp3oLP/8cjNlsF/4uXYodomBi4EAzHTqU1GjosM7DjJdGRI6V887Q2XA0Gi8pKSkhJSXF7ZoJQ1oa0atWQUAAubfaHf6OHAnkuuviSEoK4oknTnD11fa84hMnFmphqAM6dSrhtttyufXWXH79NZhNmyJ4//0IPv00gqgoG3//u30SfNgwU6XinZ5ucInCrl0h/PFHMDabwmAQeva0cO21+QwaZObss03ExjbstR7e9CjeBZYCLwBDgduAgSJyjf/Nqxm6R6FpSBw/fpycnByPcxOG1FSsrVuDUnz/fQg33hiHCCxblsnQoRXDdmjqHqsVduwIZfPmCD76KIL8/ACioqwUFQWUcTEOChL69jWTnm7g2DH7U314uI3+/c0MGmRi0CAz/fqZ/bKeIzExkbvuuotVq1Zx4YUXVvv62g49xWKPxzQGe7TZT4Dbfb0Az5doodA0FIqKivjzzz8rrJmIWrkSgoPJLzWQ/O67kdx/f0vaty9hxYp0OnXSYSgaImaz4ssvw7jzztgyc0ZOAgKEsWOLXHMMZ55p8fuajcTERObMmYPRaCQ8PJwtW7YwatSoatVRK6FojGih0DQEbDYbf/31FyJSds2EzUbc9dcjwcFk/ve/2ETx1FMxvPxyM84918iLL2b6ZfJU41u6dOmISMWJBKWEQ4fqLjJEaZFwUhOxqFU+CqVUF6XU+0qpTKVUhlJqk1Kqi9etazSnKSdPnsRkMpUVCREICCBj2TIyn32WwqIAbrwxjpdfbsa11+axcmW6FolGQtu27qdqPR33B+5EAuw92fHjx7N161aftONN0Ja1wNtAG6At8A7wlk9a12iaKBaLhbS0tDLzEmFffEHcrFmowkIIDiYlK5x//KM1X3wRxqJFWTzyyMkmEVbidGHBgmzCwsqKujPeVF0gItx2220VRMJJUVERs2bN8klb3ng9KRFZXWr/TcdCOY1G4wYRITU1lYCAgDJ5JgLy8gjIywObjZ9/Dub66+MwGgNYsSKDkSPd/7NrGi71tfjwxIkT/O9//2PdunWcOHHCY7nw8HBWrlzpkza9mcx+AsgB1mEPL341EII9lwQictInlvgQPUehqU/y8vJISko6NYFtsdiXBANYrWz5KIoFC2Jp1crKihUZdO/esMJGaBoeNpuNb775hnXr1vH5559TXFzMwIEDueaaa2jZsiU333yzX+covOlRXO14v6Hc8dnYhUPPV2g0DqxWK8ePHycsLMwe+G/fPuLmziXzP//B1K8//3mxBc88E8PAgSZeeimDli31fITGM2lpabzzzju8/fbbJCcnExMTw/Tp07n66qvp1q2bq9yKFStq7fVUGd7Eeurss9Y0miZORkYGVqvVFabD2rw5xZ06UdQsjrvuiGXTpkiuuKKAxx8/QUhIPRuraZCUlJTw1VdfsW7dOrZu3YrNZuOcc87h7rvvZuzYsYS4+eEMGzaMFStWuNZR+FIkwLuhp0exhwO3OvajgedExDezJH5ADz1p6gOj0cjBgweJjIwkMD8fW3Q0KEVmpoG5c+PYsyeEhQuzuemm3EYT1lpTdyQnJ7N+/Xreffdd0tLSiI2N5aqrruIf//gHnTp18qoOf4Xw8MbrKRD4XinVRyk1Fns48N3VtkKjacKICMePHyc4OJjA3FzaXHopzZ99lt9/D2LixDbs3x/EsmUZ3Hxz7UQiMTGR4cOHk5iY6DvjNfVGcXExH330ETNmzOD888/nxRdfpEePHixbtozvvvuOu+++22uR8CfeDD3dp5T6AtgJZAPni8iffrdMo2lEZGdnYzQaiYqKwhYcTNG4cXwVeTHXXtWGyEgbb7+dRq9ellq1Udpnfs6cOaxYsYJhdZwzW+MbDh8+7Oo9ZGVl0aZNG2699VauuuoqEhIS6tu8ClQpFEqp87GH8HgE6A28oJSaLSIp/jZOo2kMFBcXk5qaSmRJCQFZWVhbtGRxy8d4/PEYevWy8MorGcTH124RVvmFVVosGh9ms5mPP/6YdevWsWPHDgwGAxdccAHXXHMNI0aMwGComHyooeCN19MS4CoR2QeglLoCe9a5v/nTMI2msZCWlgYitL71VlTWSWb1/o517zTn4osLefrpE64EMzXF0+pbLRYNg8TERBYuXMjixYvd/h0OHDjAunXreO+998jJyaF9+/YsWLCAK6+8kvj4+HqwuPp4M5ltcE5klzrWUgcF1Gjsk4eHDh0iKiqK4i928dpiK08euIZbbsnhzjtzCPBmFtADJSUl/PLLL8ycOZO8vDyP5dq1a8c333xT84Y0Naa0iIeFhblE22g08sEHH7Bu3Tp2795NUFAQY8eO5ZprruGcc84psxDTl9RnPoquSqllQLyI9FJK9QEmAP+qtiUaTRPCarVyPCmJ5vv2cbDVcK57bCIpKYEsXZrJ5ZdXf3WuzWbjjz/+IDExke+++47vv/+eggJ7LgqlFJ4e6oYMGUJqaipt2rSp1efRVA93w4GzZs3ivPPOY+fOneTn59O5c2fuu+8+Jk2aRMuWLevZ4prjTY/ia2Ah8LLOcKfRnCIjIwNZtIhWr7zK2WG/cizkDF5+OYOzz/Yuh7WIcOjQIZcw7Nixg+xse5ygzp07c8455zBs2DCGDh3KgQMHKgw/BQcHc+aZZ7J3716UUowePZqpU6dy3nnn+e2JVWPH03Cgk3PPPZdbb72VwYMHlwkv72/qs0cRLiLfl/uwOlC+5rTGZDKRkZHBh7F387PtPIztOrPx1VQSEir/10hOTnYJQ2JiIunp6QC0bduW0aNHM2zYMIYNG1ahd+BcUOVumOPYsWOsW7eO9evX89lnn9G+fXumTJnClVdeSWxsrN++g9OVtLQ0brnlFo8iAZCUlMSQIUPq0Cr/4k2P4iNgHvCOiAxQSl0JzBGRi+vCwJqgexQafyI2G+lLn+O+PVN4fU08o0YV8dxzmURFVfxfyszMLCMMR48eBaBly5auHsM555xDhw4dvHryrGzi1GKx8Omnn7J27VoSExMJCgrioosuYsqUKQwZMqROn2ybEikpKezYsYOdO3fy/fffk5SUVGn50iJe1/irR+GNUHQBlgPnYF9HcRiYKiJ1l5mjmmih0Pian5YupdXdd5P51FO07fA34q8ax2xWoGZP4v77s3F6Nubk5LBz506XMBw8eBCA6Ohohg4d6hKGbt26+fXG/ddff7F27Vreffdd8vLy6Nq1K1OnTuWKK66gWbNmfmu3KZCcnMzOnTtd4nDs2DHA/jccPHgwQ4YMYejQoeTk5DB37twyPYv6FAmoR6EoVUkEECAi+dW2wHOdF2Ffo2EAXhWRJ8qdV47zlwBFwEwR+bGqerVQaHzJq1OXMnntXUQAhcA1EU9iNg5g5CN9mXBZOrt27XIJw2+//YaIEBYWxuDBg13CcNZZZ/nET15EXJPalf3vOs+ZTCY+/PBD3nrrLfbs2UNISAjjxo1j8uTJ9OnTp4JYVXU/8HS+Npkya3qt87qgoKAa3RiddRw7dswlCjt27CAlxb5ErHnz5gwePJihQ4cyZMgQevToUeFv6Mnrqb6od6HwNUopA3AAuBBIxh4aZLJzvYajzCXArdiFYgj2GFNVDvxpodD4itIi4aQQuClhJPvi8/n5558pKSkhODiYAQMGuIShT58+ZTPblUJEsNls2Gw213bpd0+ICEop181KKeW2V1L6mHNbKcW+fftYt24dmzdvprCwkJ49ezJ58mQmTJjgSrDk7lpvt93hycbKbK7suLtyubm5WCwWIiIiqpzEFxGSkpJcovD999+TmpoKQIsWLcr0GLp37+6VU0BV6yjqkqYoFMOwBxv8u2P/PgAR+XepMi8DX4nIW479/cBIEUmtrG4tFBpf8NPSpXS/q6xIOCkEbmjblqgJExg6dCh9+/YlJCTEJQBw6qbm/B8rfdMODAzEYDBgMBgIDAx0vZzHDAaDK/FR+Vdtyc/PZ82aNSxbtoy9e/cSFRXFtddey4033kifPn1qXX9dsnXrVmbOnMkzzzxDt27dMBgMrhDvcMqzbOfOnS5xyMjIAOzzRM7ewpAhQ/w+HFgXNEWhuBK4SESuc+xPA4aIyLxSZbYAT4jIN479L4B7RKSCCiil5gJzATp06HD2kSM1m0Kx2WzatVADQHJgIAlWz6E3jhkMmP/4g6CgINdN3nnD93STNxgMDeZmJCLs3LmTl156ifXr12MymRg2bBg33XQTV155pStUekNl69atjB8/nqKiIsLDw9mwYQM9evTgp59+Yu/evfzwww/s3LnTlQUuLi7OJQpDhw6lS5cuDeZv4Svq0z0WpdQ5QKfS5UXkjWpbUq5aN8fKq5Y3ZZz2LMc+6c7AgQNrpH4Wi4Xk5GQ6dOhAYKBXX42miVBSUoLFYsFsNlNUVERhYSF/3nknsYsXE+qmfCFw4qmn6H/GGXVtqs9QSjF06FCGDh3K0qVLWbVqFS+99BLTp09n/vz5zJw5kxtuuIHu3btXuHbr1q3MmjWLlStX+jz3gTeUFgmw54ceN24ckZGR5ObmAnZhOPfcc13i0Llz5yYnDHWFN0EBVwNdgT2A8/FKgNoKRTLQvtR+AlA+0KA3ZXyGiJCbm0tKSgrt27fXP6omiIiUEYXCwkIKCwspLranI3UOCwUGBvJX+4UE8D5j+aNMPP5C4K0pT3PdnXfWy2fwBy1atOCOO+5g/vz5fPXVV7z00ks8//zzLF26lAsuuIAbb7yRyy67jODg4DI36fHjx/ssm5rVaiUnJ4eTJ0+6XtnZ2RX2Dxw4wPfff+8a4it9fUFBAQsWLOD666+nefPmZGRkEBAQUGY4SlN9vHGP/R04S3w8RqWUCsQ+mT0aOI59MnuKiPxWqsw47Gs4nJPZz4vI4Krqrukchdls5uDBg9hsNhISEmjRokW169A0HESE4uJiLBYLJpPJJQpWx3CSUxSCgoIqeLNs2RLObbe1JCpsCH8v3sfKYqPL6+mtKU9z3ZqmIxKeSEtLY+XKlbz88sscOXKE+Ph4xowZw4YNGyrNz2w0Giu90Xvad/YEPBEVFUWLFi1ISUlxCbs7Onbs6FrrYDabSUtLIy8vj9DQ0Bp7RzUW6nMdxTvAbVVNINcEh1fTs9jdY18TkceUUjcCiMhLDvfYF4CLsLvHznI3P1Ge2gpFWFgYRUVFdOvWjdBQdwMPmoaGiGCxWLBYLBiNRgoLCzEajS5RCAgIcIlCVXNQ21ZnEbroKf7Vfhi/HL2DJ554gjEGA/H33kvmU0/Rvwn1JLzBarXyySef8K9//ctjwiSlFDExMRQVFWEymTzWZTAYaNGiBS1atCAmJsa1XdV+8+bNXTe/8sNOpXGXL1pEKCws5Pjx4157RzVW6lwolFLvYx9iigL6Ad8DriA2IjKh2pbUEbUVisjISMxmMwEBAXTp0qVBx4k/HbHZbC5RcM4nGI1Gl/uoUso1wVzdG8JHH4XzybwdrAiYw6XhxRi7dObLL78kJibGT5+m8dCpUycqcxKJjIzkpptuqvTmHxkZ6ZMhIHdi4U4kSmOz2cjKyiI9Pb3JDkfVx2T2kmq31IQICQmhoKCAtLQ02rVrV9/mnNaYzWaMRiNFRUUUFRWVGfYwGAwEBQURERFR63/6jz8K47bbWtG3/1juOeMidry9hg2LFtG8efNafoKmwcqVKyt9kt+8eXOdTWyPGjWKLVu2lPF6qmquJCAggFatWtGsWTNSU1NPm+EoX+DxcUtEvhaRr4FLnNulj9WdifVHREQEJ06cqHLsVOM/8vPz+fPPPzl27Jjr7xAZGUlUVBRRUVGEh4cTFBRUa5H4Ygt0umUmczp9zP33b+P1d9/iyiuvZMyYMU3uqbOmOG/O4eHhZY57c5P2pz0dO3asVvvBwcF06NCBzp07Y7PZyM/PrzAxrimLN/3yC90ca7ABAX2JUoqIiAiOHTuG2exd6GiN78jNzSUpKYng4GCioqIICwsjMDDQ5zfuTz8N48H5oXQOSeHeGw6zePEiIiMjuffee10rljV2yotFfYlEaXuSkpKq3b5SisjISLp160br1q1dvdX6WlfW0PEoFEqpm5RSvwA9lFJ7S70OA3vrzsT6xbmQKjk5WT911CHZ2dkcOXLE1WPwF59/GsK8W1rRqlcM5m83sCXUxo4dO5g3bx5nnnmm7k24oaZP8g0R53BU9+7diYiIID8/v1KPqtOVyiazmwExwL+Be0udyheRk3VgW43xxWR2efLz84mLi2s0OW4bKyJCVlYWKSkpRERE+NWR4MsvQrHNfZDo5vC3L/+PwKBCxowZQ/Pmzfnwww/p1KmT39rWNEwKCgoatXeUvyazK/sWDEAecAuQX+qFUuq0W2AQERFBeno6+fk+C56rKYeIkJGRQUpKCpGRkX4Via++CuOmm+KgZTNGTAwiupnw4osvkpqayr333kvbtm391ram4eIcjmrTpk2jGI6y2WyYzWYKCgooKCjwWw+4Mq+n3ZwKl1G+dQG6+MWiBorTne7YsWN069ZNe0r4GBEhNTWVEydOEBkZ6dcnua+/DuP2uVF071HMoDdvxRhtJSkpiVdffZVLL72UsWPHeoz8qmn6BAQEEBsbS3R0NGlpaeTk5BAaGlrvvwmbzUZxcXGZobGAgAAiIiJo2bIloaGhhISE+OXe5FEoRKSzz1tr5AQFBVFcXMzx48fp2LGjHr/2ETabjePHj5OdnU1UVJRfv9dt20L5ds4Wfgl4nPyn3yKiWRygePTRRwkKCuLOO+/U6UM1wCnvqBYtWnD8+HHy8/PrbDjKW1Hwh3OHO7wNChgDdINT8dFEZJu/jGrIhIeHk5eXR1ZWlr6h+ACr1cqxY8fIz8/3u0hs3x7K3LlxTGzfg8gze2Pr1ByAL774gi+//JIFCxbQu3dvHRBSUwbncNTJkydJS0vz+WK98qIgIhgMhnoTBXd4ExTwOuB27AH59gBDgUTgAr9a1oCJjIwkJSWF8PDwCj7lGu8pKSnh6NGjFBUVERUV5de2vv02lHuvC6BzlxLuXtOBvBbPAnYHhkcffZQuXbowbdo0vQJb4xZfDUc1BlFwhzePTrcDg4AdIjJKKfU34GH/mtWwCQgIIDQ0lGPHjtG1a1f9BFoDiouLSUpKori42K2XmS/57rtQ/jP7OPtLRpEx60kCWpxaBvTqq69y5MgRXn75ZTp27NjovFw0dUv54aiCggLCw8Pd/m4aqyi4w5s7nElETI4YOiEi8odSqoffLWvgBAcHU1hYSGpqKgkJCQ3+D92QMJvNJCUlYbPZ/N4j27EjlDlz4ujRMRRjn8sIHDUQ52qYlJQUXnzxRcaOHcsFF1xAdHS0X23RNB3KD0c5IxE3BVFwhzdCkayUag5sBD5TSmXjx5wQjYnw8HCys7OJiIjQIcm9xGQycfjwYZRSfs+gtmNHCA/MstClXSGvrs2jMPbRMucff/xxbDYb8+fPp02bNo3yH1hTf5QejsrIyMBqtTYJUXBHlUIhIpc7Nh9SSm0FmgEf+9WqRoIzxMfx48cJDw/XIcmroKioiMOHDxMUFOR3V8OdO0O4eVY0e0p6E9G1N/mxz5U5n5iYyAcffMC8efM466yzdKgOTY0JDg4mISGhvs3wK956PQ0HuonISqVUK6AdcNivljUSnNFLjx07pkOSV0J+fj5JSUl1Eq3zhx9CmD07njbtSiiZOg/zgK5lzhcXF/PQQw+RkJDAtGnT9Gp7jaYKvPF6WgQMBHoAK4Eg4E3gXP+a1ngIDQ2loKCA9PR0vaLXDbm5ua64Tf6e+N+1K4RbZ4Zzbotf+dfaGALirsBSrsybb77JgQMHeO6552jXrp3uCWo0VeCNi8flwATsWSARkRTsyYw0pdAhyd2TlZXFkSNHiIiI8LtI7N4dwsyZ8SxTN/NB4QXEh+VUKJOZmckzzzzDeeedx6hRo4iLi/OrTRpNU8Cb/1yLiIhSSgCUUnow1w1KKcLDw0lOTiYsLKzel/vXNyJCZmYmaWlpfg/JAfDTT3aRaNXKSqelt5Cdei7iZm3G4sWLMZlMLFy4kPj4eB2KRaPxAm/+e99WSr0MNFdKXQ98DrziX7MaJ87Um6d7SHIRIS0trQ5FIpiZ01txVfj7vLU2lZj+bSm6pGJurZ9++ol33nmHmTNn0rVrV1q2bOlXuzSapkKV/8EisgR4F/gf9nmKB0XkP/42rLESFhZGYWEhmZmZ9W1KvWCz2UhJSSEzM5OoqCi/i8SePcHMmNGamSFreC1jIp3+/NqjXQ899BBxcXHMnj2bNm3aaMcDjcZLvBo0FpHPgM/8bEuTwRmSPCIiwu+rjhsSNpuN5ORkcnNz/R63CWDv3mCmT29N8+ZWrl47koy9L2IaPtxt2bfffpu9e/eyZMkSYmJiaNasmV9t02iaElU+7imlrlBKHVRK5Sql8pRS+UqpvLowrrFSOiT56ZIty2q1cuTIEfLy8oiMjPS7SPzySzDTpsUzJWg9by/7nbYJYh9uctNubm4uixcvZuDAgVx44YW0adNGh+rQaKqBN/8tTwETRKSZiESLSJSI1CrWgVKqhVLqM4cAfeaITlu+THul1Fal1O9Kqd+UUrfXps26JigoyDUM05ATn/iC4uJiDh8+TFFRUZ2IxK+/2kWia3gKLxbO5swNz1VafunSpeTk5PB///d/REZG+j0AoUbT1PBGKNJF5Hcft3sv8IWIdAO+oGyqVSclwF0icib2iLW3KKXO8rEdfiUiIoLc3FyysrLq2xS/YbFYSEpKwmw218nq5t9+C+baa+OJiBCee1uRvu4tcu6+22P533//nTfffJMpU6bQpUsXWrdu3WTCKmg0dYXHOQql1BWOzV1KqfXYYz2ZnedFZEMt2r0MGOnYXgV8BdxTuoCIpAKpju18pdTv2FeE76tFu3VOZGQkqampTTIkudls5vDhw4hInYjEvn1BXHttPBcEbuPh25OJan8elvb9PJYXERYtWkSzZs246aabiImJaXJ/A42mLqhsMvvSUttFwNhS+wLURijiHUKAiKQqpSpd9aSU6gT0B3ZWUmYuMBegQ4cOtTDNtwQEBBASEtLkQpIbjUYOHz7smo/xBxs3RrBkSQwpKQZatbJSUBBAs2grK9r8k4jXskiZ9AFU4rn0/vvv88MPP/Cvf/2LqKgovbhOo6khlaVCnVWbipVSnwOt3Zz6v2rWE4ndNXe+iHicRBeR5cBygIEDBzaoSYHg4GAKCgqaTEjywsJCDh8+THBwsN8WFm7cGMH997dkiPFrvmEWszJW8hUjueOOPPL+sYyCoqJKRaKwsJDHH3+c3r17c+mllxIbG0tISIhfbNVomjp+e7wVkTGeziml0pVSbRy9iTZAhodyQdhFYk0th7rqnYiICLKzs4mMjGzUWdTy8vI4evSo35K4O1myJIYhxq/ZwngiKOIDxrGZCdyz8nWuuy4aaxW5I1544QXS09N54YUXXOGgNRpNzagvH8HNwAzH9gxgU/kCyv7YvQL4XUSW1qFtfqF0SHKTyVTf5tSInJycOosA2/34NpdIAIRj5GrWc37K21Vee+jQIVasWMGkSZPo0aMHrVu31qE6NJpaUF9C8QRwoVLqIHChYx+lVFul1IeOMucC04ALlFJ7HK+KcRkaEQaDgcDAQI4dO4bVaq1vc6rFiRMnOHr0aJ0E9wtNTGQLl7pEwokCXlY3EZqY6PFaEeGRRx4hNDSUu+66i8DAwEbdg9NoGgLehBkPASYBnUqXF5FHatqoiGQBo90cTwEucWx/g/3e0KRobCHJ6zq4H0DU7XcTXk4knIRLEcELF5L8zTduz3/++ed8/fXX/POf/yQyMlKH6tBofIA3//WbsLuzlmAPNe58aWqIMyR5Xl7DXuBe18H9wJ6Z7qr81ylS7t1YbWFhnFi82O05k8nEo48+Srdu3bjmmmsIDQ3VoTo0Gh/gzRhCgohc5HdLTiOcIcmPHTtGt27dGlRIcqvVislkorCwkJycHEwmU53EbQL49ttQrr8+jrbtWnLkztfosWAWAUaj67wtLIyMFSswDRvm9vrly5dz7Ngx3nzzTaxWq86DrdH4CG8eEb9TSvX2uyWnGQ0lJLnNZqOoqIgTJ05w6NAhfv/9dw4dOkRmZiYBAQFER0fXyc3266/DmDMnjg4dSnhrbSotIkxkvPoqNscajapEIjk5mWXLlnHxxRdz9tlnExUVpfNgazQ+whuhGA7sVkrtV0rtVUr9opTa62/DTgecIclPnDhRZ22KCCaTiezsbI4cOcLvv//OX3/9RVpaGiUlJURERLhusnW1OPDzz8OYOzeOrl2LWbs2jQ57PiF+5kwCCgrIWLGCknbtKhUJgMceewyA+++/H4vFokN1aDQ+xJs7wcV+t+I0JiIigrS0NCIiIvz2BGyxWDCZTOTl5ZGXl+fyuAoODiYsLKxeI6l+9FE4t93WirPOsrBqVTrNm9swjhlD5vPPUzRmDAQEeJy4dvLtt9/y8ccfc+eddxITE0Pz5s39tlpcozkdqSzWU7RjJXR+Hdpz2uEMgXH06FHOOOMMn/j7l5SUYDKZKCgoIDc3F4vFglIKg8FAaGhogwmxvXlzBHfeGUu/fmZeey2d6EgrqqAIiYyk8NJLq64Ae+Tahx56iA4dOnDddddhtVpp1aqVny3XaE4vKutRrAXGA7uxx3Yq3Y8XoIsf7TqtCAoKwmKxkJKSQocOHao9ZGKz2VwT0Lm5uZhMJkQEg8FAcHBwgwyr/b//RXD33bEMHGhmxYp0IiOFqJWraPbqq6Ru2IA1Pt6relatWsWff/7JK6+84hKJhuQcoNE0BSqL9TTe8d657sw5fSkdkryqcBMigtlspqioiNzcXAoLCxERAgICCA4OJiIiokGPz69bF8n997fknHNMLF+eQXi4PTSXpV8/isaMwepl8L7MzEyee+45Ro4cyciRI7FYLDoPtkbjB5pGKNMmgjMkeURERJkxdhHBYrFgNBrJz88nPz8fq9WKUqpRCENp3ngjikWLWjJiRBEvvZRJaOip+I3m/v0x9+/vdV1PPPEEZrOZBx54AJPJRJs2bZpMdF6NpiGh/6scrFkD998fzLFjPWnb1sqCBdlMnFi36wqdIcmPHj1Kp06dsFgs5Ofnk5eX50qpGhQU1KDmGarDq69G89hjLRgzpogXXsjAGcw1+pVXUCYTubfcAl5+rt27d7NhwwZuuukm2rdvj81m06E6NBo/oYUCu0jMnQtFRfan8uPHA7n/fvsQRl2LRXBwMIWFhRw8eBDANc8QGhpap3b4mv/+txmLF8dw8cWFPPdcJq45exGC9+9HFRW5zXftDqvVykMPPUTr1q255ZZbMBqNdOzYsVGKp0bTGPBKKJRSBiCesrGejvrLqLrm//4PisqFFjIaA1iyJKbOhQJoUgvFROD555vx7LMxXHZZAUuWnKDM6JBSnFiyBMxmr4Vi/fr1/Prrrzz//PMEBgYSEhJCdBVhxzUaTc2p8hFMKXUrkA58BnzgeG3xs111ylEPkpeSooPJ1QYRWLKkOc8+G8OVV+bz9NNlRSJi0yYMaWn2HS+TCmVnZ7N48WKGDBnCuHHjMJvNOlSHRuNnvOmr3w70EJGeItLb8erjb8PqEk+ZU9u0aVyhwBsSIvDYYzH897/NmTw5nyefzCqTkC4gN5eWDzxAsxdeqFa9S5cuJS8vj0WLFmE2m2nWrFmT6oFpNA0Rb4TiGJDrb0Pqk8ceg/AKwUqFyEgbRqN+Uq0uNhssWtSCFSuaMWNGHo89llVhjtrWrBkpmzaRc/fdXte7b98+1q5dy7Rp0/jb3/5GSUkJ8V6ut9BoNDXHG6E4BHyllLpPKXWn8+Vvw+qSqVNh+XLo0EFQSmjXroTJk/M5eDCI666Lc01ya6rGZoP772/J6tXRXH99LosWnaww9RB4+DAAJZ07Y/NybkFEWLRoEc2bN+eOO+7AaDTSokWLRj/Jr9E0BrwRiqPY5yeCgahSrybF1Klw4ICFvXt/45tvknn88ZM8/fQJduwIZfbseAoLtVhUhdUKd98dy/r1UdxySw733ZddQSRCfvyRdmPGEPH++17VmZiYyPDhw1m8eDG7du3i7rvvJioqCpvNRpyXC/M0Gk3tqNLrSUQeBlBKRdl3pcDvVjUQLr+8EIMB7rgjlpkz41m50h5qQlORkhK4665YNm+O5I47srntNvejlZazziJn/nyKLrigyjoTExOZM2cORqORZcuW0aVLF6666iqKioqIi4vTebA1mjrCG6+nXkqpn4Bfgd+UUruVUj39b1rDYMIEu9//Tz+FMGNGPPn5umdRHosFbr21FZs3R3LPPSc9igQ2GxIaSu6ttyJVTECXFgknx48f57vvvsNgMOhQHRpNHeLN0NNy4E4R6SgiHYG7gFf8a1bDYvz4Iv7zn0z27g1h+vTW5OVpsXBiNsPNN8fx8ccRPPDASW680X1619Bt22gzcSKGlJQq63QnEva2zFx//fUcPHhQ58HWaOoQb4QiQkS2OndE5CvgtPNHvPjiIl58MYPffgtm2rTW5ObqVcAmk+KGG+L44otwHn00i9mzPecAVzYbtogIbF70BBYuXFhBJE61aeLWW2+tsc0ajab6eOX1pJR6QCnVyfH6J3DY34Y1RMaONbJsWQZ//BHMtdfGk5Nz+opFUZFizpw4tm0L44knTnDttZWnLTGOHEn62rWIFwvrFi9e7DFUeFhYGK+//npNTNZoNDXEmzvdbKAVsAF4z7E9qzaNKqVaKKU+U0oddLx7jOamlDIopX5SSjWI1eCjRxt56aUMDhwIZurUeE6ePP3EoqBAMXNmPDt2hLJkyQmuvtqzf0PYF18Q+e679hV4Xq6ezszMdEXHLU1oaChbtmxh1KhRtbJfo9FUjyrvciKSLSK3icgAEekvIreLSHYt270X+EJEugFfOPY9cTvwey3b8ymjRhlZvjydv/4KYsqU1mRlnT5ikZenmDEjnh9/DOG55zK54orKY2FFbthA1MqVdreoKhARli9fzu23387AgQNZvny5K9x6aGgoGzZs4AIvvKU0Go1v8XiHU0o963h/Xym1ufyrlu1eBqxybK8CJnqwIQEYB7xay/Z8zogRJl59NYMjRwKZMqU1mZlNXyxycgKYNq01v/wSwosvZjJ+fFGV12T+5z+kv/46VOHKarVaeeSRR/j3v//NuHHjWLVqFWPGjGHFihW0bduW119/nYsv1unbNZr6oLJ1FKsd70v80G68iKQCiEiqUsrTyqlngbtpoAv8hg838dprGcyZE8fkya1ZuzaduLimGR/q5MkApk2L588/g1m2LIPRo91PNjsJ+eknLGecgURFYasih7XZbGb+/Pl8/PHHzJkzh/vvv98VMvzss8/m008/5YwzzvDZZ9FoNNXD42OwiOx2bPYTka9Lv4B+VVWslPpcKfWrm9dl3himlBoPZJSyo6ryc5VSu5RSuzIzM725xCcMG2bitdfSSU0N5JprWpOW1vTcNjMzA5g8uTV//RXEK6+kVykSymgk7rrriL3vvirrzsnJYdq0aXz88cf885//5J///KdLJCwWCxaLhS5duhDiZXRZjUbje5RI5SuNlVI/isiAcsd+EhHvc1ZWrHM/MNLRm2gDfCUiPcqV+TcwDSgBQoFoYIOIXFtV/QMHDpRdu3ZV2y6z2czBgweJjIys9rU//BDCrFnxtGplZe3atEYdeXbjxgiWLIkhJcVAfLwVqxUKCgJ49dUMzjnH5FUdIT/+iDUmhpLOnlOuHz9+nJkzZ3L06FGefvppxo8f7zpXXFyMxWKhc+fOhFeM2KjRaHyMUmq3iAx0d66yOYrJSqn3gc7l5ie2Alm1tGkzMMOxPQPYVL6AiNwnIgki0gm4BvjSG5GoLwYNMvPGG+mcOGHgmmtac/x44+xZbNwYwf33t+T48UBEFGlpgWRmGpgzJ9crkVAFdg8o84ABlYrEvn37mDRpEunp6axataqCSJjNZi0SGk0DobIZ2O+Ap4E/HO/O113ARbVs9wngQqXUQeBCxz5KqbZKqQ9rWXe9MWCAmTfeSCM72y4WycmNL9PskiUxGI3lfxaK996repooMDmZhPPPJ2Jz5b4O3377LVdffTVKKd555x2GDh3qOldcXIzJZNIiodE0ICqbozgiIl+JyLBycxQ/ikjVvo6VICJZIjJaRLo53k86jqeIyCVuyn8lIuMr1tTw6N/fwptvppGXF8DVV7fm6NHGIxZFRcpjT8ibbH+2iAiKRo/GPGCAxzIbN25k1qxZtGvXjg0bNtCjx6kRR6dIdOnSRScj0mgaEN4EBRyqlPpBKVWglLIopaxKKc+xGjT06WNhzZp0iooUkye3JimpYYtFaqqBJ5+M4ZxzEgD3i+Latq16zsUWE0PW4sWUJCRUOCcivPTSS9xxxx2cffbZvP3227Rp08Z1vqSkBLPZTKdOnbRIaDQNDG+c/18AJgMHgTDgOuA//jSqKdCrl4U1a9IwGu1icehQwxOLn34K5rbbYjnvvASWL4/mnHNM3HZbNmFhtjLlwsJsLFjgeY1l4KFDtLrhBgwZGW7PW61WHnroIZ588knGjx/P66+/TnSphEUlJSUYjUY6duxIVFSD9ITWaE5rvLp7icifSimDiFiBlUqp7/xsV5PgrLOKWbs2jWuvbe1aZ9G1a3G92lRSAh9/HM5rr0Xz00+hREXZmD07j+nT80hIsPcaOncucXk9tW1rZcGCbCZO9LwCO3j/fkJ++cUepqMcJpOJ+fPn88knn3Dddddx3333udxf7fZokdBoGjreuMduA8ZgXx2dBqQCM0Wkr//Nqxn14R5bGQcP2kN9AKxdm0a3bnUvFrm5Abz1ViSrV0eTkhJIp07FzJyZx6RJBb5JxmQ2Q7m1Djk5OVx//fXs3r2bBx54gFmzyoYIs1qtFBYW0qlTpzI9DI1GU/fUyD22FNMAAzAPKATaA5N8Z17Tp1u3YtatS0MpYfLk1uzfX3eZ2f76K5AHHmjBsGEJPPlkCzp1KubVV9P54ovjzJiRXyuRCNq/n9Bt2+w75UQiOTmZK6+8kr179/Kf//zHo0h07NhRi4RG08DxJhXqEcemEXjYv+Y0Xbp2Leatt9KYOrU1U6a0ZvXqNM46yz89CxH45ptQVq6MZuvWcIKDhYkTC5g5M48zz/Rdm81feIHQHTtI/vprpJQr6759+5g5cyZms5k33niDIUOGlLmutEg0a9bMZ/ZoNBr/4FEolFJvi8g/lFK/ABUeO0Wkj18ta4J07VrCunVpTJnSmqlTW/Pmm+n07GnxWf0mk+K99yJ4/fVoDhwIJjbWyh13ZDNlSj6xsbaqK6gmJ5YsIfDQoTIisX37dm6++Waio6N588036d69e5lrnCLRoUMHLRIaTSOhsh7F7Y73RrF+obHQqZNdLCZPbs2UKfGsXp1Onz61E4v0dAOrV0exdm0U2dkGzjrLzJIlmYwfX1h+RMgnGI4fxxoXh4SEUHzmma7jGzZs4J577qFr1668/vrrtG7dusx1TpFo3749zZs3971hGo3GL1S24C7VsXkFUOJYgOd61Y15TZMOHUpYvz6N6Ggb117bmj173Gdzq4q9e4OZPz+W4cMT+O9/mzFokIl161LZsiWVSZP8IxIUFxM/cyZxN9/sOiQiLFu2jLvuuotBgwbx9ttvVxAJm81GYWEhCQkJxMR4zFOl0WgaIN64x0YDnyqlTgLrgHdFJN2/ZjV9EhJKXHMW06e35vXX0xkwwFzldSUl8NlndvfWXbtCiYy0MX16HjNm5NOhQ60WzHtHUBA5CxZgc7iyOtdIvPnmm0yYMIGnnnqqQqRXm81GQUEB7dq1o0WLFv63UaPR+JQq3WNdBZXqA1yN3eMpWUTG+NOw2tDQ3GMrIyXFwNSprcnMNPD66+kMHOheLPLyAli/PpJVq6I5fjyQ9u2LmTkzn6uuyicqygfurd5gs0GpNRAmk4nbb7+dTz/9lLlz53LPPfeUWSNhv+SUSLRs2bJu7NRoNNWmMvfY6iwXzsC+jiIL8JRoqFGjvvqK7rNmcfLppzENG1YnbbZta+Wtt+wT3DNmxDN7di7vvRflWuw2Y0YuyclBvPtuJEVFAQwZYuLBB08yenQRhjoMUKsKC2l99dXkzptH0UUXkZ2dzfXXX8+PP/7Igw8+WMH9FU6JRNu2bbVIaDSNGG8W3N2EvSfRCngXWC8i++rAthpTox7F1q3I+PGooiJsYWFkrFhRZ2IBkJFh4NJL25CRYaBsvCUhIAAmTixk1qw8evXynZdUdTBkZBA7fz458+fzV9u2zJgxg+TkZJ599lm3KUpLi0RsbGw9WKzRaKpDbXsUHYH5IrLHp1Y1JLZuBYdIAAQYjcTNmVOnYhEXZ3WM6pQPyqeIiyvh6adP1IkdnrDGxZG+Zg2/7dvHrCuuwGKx8OabbzJo0KAKZZ0i0aZNGy0SGk0ToMqV2SJyLxCplJoFoJRqpZTynJGmseEQCRwi4cQpFqGJiXVmSnq6+7EkT8f9SWhiIgnDhxP2+efEPPYYKj+fbdu3c/XVVxMUFMS7777rUSTy8/Np3bq1FgmNpolQZY9CKbUIGAj0AFYCQcCbwLn+Na2OmDWrgkg4CTAaiV24kORvvqkTU9q2tXL8eMU/iTchvn1JaGIisTNnEmixEHvTTQQAmyMjmfPCC3Tr1o2VK1cSHx9f4brSPYlWrVqhlPuQ5RqNpnHhTayny4EJ2OM8ISIpQNMJ87lyJXjIpGYLCeHEU0/VmSkLFlQ/xLevKS0SAIaSEiwivPbsswwZMoT169e7FQkRoaCggPj4eC0SGk0TwxuhsIh9xlsAlFJNK6vMqFGwZUsFsbCFhKCsVoL37q0zUyZOLOTxx7No164EpYR27Up4/PGsSkN8+5LyIuEkxGrlw4AA3po7120ocBEhPz+f+Ph44uLitEhoNE0Mb7yeFgDdsOe2/jcwG1grIg02eZFPvJ5eeYWgAwcovPRSbLGx9pVugQ0v+ZAvaTVoEBEnPE+aF8bGkvnDD2WOOUWiVatWtG7dWouERtNIqVWYcRFZgt0t9n/Y5ykebMgiUWNGjaJ4wwYsbdrYvZ3OPZf8WbPsIiFC3M030+Lhph0893mLBfezNfZxx/IrJbRIaDSnB95muPsM+MzPttQ7MnIkBz/9tOLKbKuV4s6dsboZm28qZO3cyd15ebwB/AMoPb5YCEwKDmbG88+7jjnnJGJjY7VIaDRNnMrCjOfjJry4ExE5fbLNBAaSfd99rt2QHTsI++47cm69FYLqLgmRr1EmE9nvv8+/d+1iw4YNjDYYKDn3XNYlJrKhuJgISonE668zzLGmxNmTaNmyJW3atNEiodE0cTwKhYhEASilHsEeumM19tVgU2lKXk81IGzbNiI++ojcG29EGqlQ/PHHHxTNm8f4v/7ih+Bgpk6dyty5c2nbti2JiYlMmjmTlywWbnQjEgUFBbRs2ZK2bdtqkdBoTgO8mczeKSJDqjpWrUaVagGsBzoBScA/RKSCD6hSqjn2XN29sPduZotIlSvg6iIoYEBeHrboaCgpIfzjjykaNw4awU3zwEcfsWrNGtZ++y0dw8O57YILOO/BB2nVqlWZcomJiSxcuJDFixdXEImYmBjatWunRUKjaULUNoSHVSk1FXuIcQEmA7VdAXYv8IWIPKGUutexf4+bcs8BH4vIlUqpYMD9god6wObI8xzx/vu0uvNO0po3xzR8eD1b5R4RYefOnbzy/POsS0ykMCiIuPnzmTFjhscEQsOGDeObUgsNnSLRvHlz3ZPQaE4zvBGKKdhv2M9hF4pvHcdqw2XASMf2KuArygmFUioaOB+YCSAiFqB+IuJVQuHEiVhbtHCJhCE9vcFMeosIX3/9NVufeoo3fv+d2NhYPrjqKobeeCP9u3SpVl0FBQU0a9aMdu3aVQglrtFomjZVCoWIJGG/sfuSeGcGPRFJVUq5C1veBcgEViql+gK7gdtFxO3qM6XUXGAuQIcOHXxsbiUohWnECMAuEm3//nfyrr+e3FtuqTsbymGz2fj000958cUX6f3rr6wFRs2YwdB77yU0NLTa9eXn59OsWTMSEhK0SGg0pyF+W0GmlPocaO3m1P95WUUgMAC4VUR2KqWewz5E9YC7wiKyHFgO9jmK6ltce6wtWpA3ezaFzrDbInU6b1FSUsL777/Pay++SOFffyGdOnH2v/5FenExI6dOrZGHVkFBAdHR0VokNJrTGL8JRWUZ8JRS6UqpNo7eRBvsSZHKk4w9k95Ox/672IXCr9hsNkSkZmPwQUHk3naba7fFQw8hoaFk33uvXwXDbDazYcMGXnrpJY4ePcrOsDA6t21L7iefEBgcjLEGdToD/DVr1oz27dtrkdBoTmPqKybFZmAG8ITjfVP5AiKSppQ6ppTqISL7gdGAXxMmBQcH06pVK06cOEFQUFCNhmlc2Gz2l9XqN5EwGo2sW7eO5cuXE5KWRkzv3vzfyy/TPjQUa3AwgcHBNa63pKSENm3a0LJlSy0SGs1pjtdCoZQaCjwOhACLRWRjLdp9AnhbKTUHOApc5WijLfCqiFziKHcrsMbh8XSIilEkfIpSijZt2tC8eXOOHz9OXl4eERERGGqSczQggJOPPmoffgKCDh4kbOtW8ubMobY5TPPz81m9ejWvvfYaWVlZTO7Vi9VZWWRPnkzB2LG4z7pdNcXFxZhMJqKjo2ndujUhISG1slOj0TQNKluZ3VpE0koduhN7uHEFfAdsrGmjIpKFvYdQ/ngKcEmp/T3Yc2HUKWFhYXTt2pWcnBxSU1MREcLDw2s2HOW4JuK994hav56CSZOw1TB/dHZ2NitXruT1118nPz+fy4cOZfKyZQwaOJC8//wH4wUX1Khem81GUVERBoOBDh06EB0drd1fNRqNi8rGFF5SSj2glHKOv+Rgd4u9Gsjzt2H1jVKKmJgYunfvTvPmzcnPz8dsrumzOuQsXEjK5s12kRCxZ87zsNgxMTGR4cOHk+jIrpeRkcFjjz3G8OHD+c9//sM555zDgUmTeOfPPxncowcoRe5tt9XILddkMrliNnXv3p1mzZppkdBoNGWoLITHRKXUpcAWpdQqYD52oQgHJtaJdQ2AwMBA2rVrR0xMDCkpKeTn5xMeHl794SilsLZrB0DYF18Qf/31ZLz8MkVjx5YplpiYyJw5czAajcyaNYvzzjuPbdu2UVJSwqRLLmHu3Lmc0bs3wXv3ktujBxIWVqPPZbVaKSwsJDIykk6dOtVuPkaj0TRpvAnhYQBuBsYBj4nI9rowrDbUNIRHVdhsNrKzs0lLs4/I1Xg4ymolYtMmCi+7DAwGVH4+EhVVRiRKM3LkSB5ZuJChN95I4aWXkrNwYY0/g4hQVFSEUoq2bdvqHoRGowFqmI9CKTVBKfUN8CXwK3ANcLlS6i2lVFf/mNqwCQgIoGXLlmWGoyyWGiwWNxgovOIKl0i0HTeODEdIDaPRyEjgMPal6xHAzp07Sc7NpWDSJEznnVdj+81mM/n5+TRv3tz1GbRIaDSaqvDYo1BK7QWGAWHAhyIy2HG8G/CoiFxTZ1ZWE3/1KMpTWFhISkoKJpOpRsNRhYWFfPnhh7R+/nmeSU4mEbs4bMEuECagGOgDFLdrVyb2UnWwWq0UFRURGhpKu3btCPeQI1yj0Zy+1DQoYC72XkQYpRbEichBx/HTnoiICLp27Up2djapqakopaocjrJYLGzbto3Nmzfz2WefYTKZaNu2LWdfeiljP/6Y94uLca5+CAUMwNCgICYsXlxt+5zDTCJC27ZtiYmJ0WsiNBpNtalMKC7HHim2mNoHAWyyOIejoqKiSE9PJzs7m9DQUIJLLXaz2Wx8//33bN68mQ8//JDc3FxiYmK48sormTBhAmeffTbhO3cS+8knFf4gQcBqpTiBvYfhLcXFxRiNRmJiYoiPjy9jj0aj0VSHKiezGyN1NfTkjoKCAlJSUjCbzSQlJfH+++/z/vvvk5aWRnh4OBdeeCGXXXYZw4cPJ6hU7KWE4cMJPH7cY70l7dqR7MXQk81mo7CwkODgYNq1a+dVbg2NRqOpbT4KTTVIS0tj3bp1vPnmmxw8eJDAwEBGjBjB/fffz+jRoz3OD5xYvJi4OXMIMFaMzGQLC+OEF0NPztAb8fHxxMbG6mEmjUbjE7RQ+IC0tDTWr1/P2rVr+f777wEYMWIE8+fP59xzz0UpRVhYWJkeRHlMw4aRsWJFBbGwhYWRsWIFJkeWOXc4h5mio6Np06aNDr2h0Wh8ihaKGpKbm8uGDRtYu3YtX375JTabjf79+7N48WKuvvpq2rdvD9gnlAsLCzl+/Dj5+flERER4fNIvLxZViUTp0BsdO3bUoTc0Go1f0HMU1cBkMvHBBx+wdu1aPvjgA8xmM127dmXKlClMnjyZM8880+O1NpuNrKws0tLSCAwMJDQ01ONNPTQxkdiFCzmxeLFHkTCZTFgsFlq1akVcXFzNAhdqNBqNg8rmKLRQlGLr1q3MmjWLlStXMmrUKMCeDGjr1q2sXbuWDRs2kJeXR3x8PNdccw1Tpkxh0KBB1XqKN5vNpKWlkZubW+VwlDtKh95o06YNYTUM4aHRaDSl0ZPZXrB161bGjx9PUVER48ePZ/Hixezfv5/169eTnp5OdHQ0kyZNYsqUKYwcOZLAwJp9dSEhIXTo0IH8/HyXd1R4eHiVE8/ONREA7du316uqNRpNnaGFgrIiAVBUVMQtt9xCUFAQl156KVOmTGHcuHE+C5ynlCI6OpqIiAiysrJIT08nMDDQY+/AYrFgMplo2bIlcXFx1e6FaDQaTW047YWivEiUJjAwkHnz5rmGoXyNwWAgLi6OZs2akZKSQl5eXpnhqNKhN8444wwdekOj0dQLp/0cRadOnThy5IjH8x07diQpKclHlnlGRMjLyyMlJYWSkhICAgIQEVq3bk2LFi30mgiNRuNXahQ99nRh5cqVHp/Uw8PDWblyZZ3YoZSiWbNmdO/e3dXL6N69u144p9Fo6p3T/g40atQotmzZUkEswsPD2bJli9+GnTxhMBiIj48nISFBx2fSaDQNgtNeKKCiWNSXSGg0Gk1DRAuFA6dYdOzYUYuERqPRlOK093oqzahRo+pk4lqj0WgaE/XSo1BKtVBKfaaUOuh4j/FQ7g6l1G9KqV8dKVh9s5BBo9FoNF5TX0NP9wJfiEg34AvHfhmUUu2A24CBItILe7I3nVlPo9Fo6pj6EorLgFWO7VXARA/lAoEwpVQgEA6k+N80jUaj0ZSmvoQiXkRSARzvceULiMhxYAlwFEgFckXkU08VKqXmKqV2KaV2ZWZm+slsjUajOf3wm1AopT53zC2Uf13m5fUx2HsenYG2QIRS6lpP5UVkuYgMFJGBrVq18s2H0Gg0Go3/vJ5EZIync0qpdKVUGxFJVUq1ATLcFBsDHBaRTMc1G4BzgDeranv37t0nlFKe43I0DmKBE/VtRANBfxdl0d9HWfT3cYrafBcdPZ2oL/fYzcAM4AnH+yY3ZY4CQ5VS4YARGA14FcBJRBp9l0IptctT3JXTDf1dlEV/H2XR38cp/PVd1NccxRPAhUqpg8CFjn2UUm2VUh8CiMhO4F3gR+AXh63L68dcjUajOX2plx6FiGRh7yGUP54CXFJqfxGwqA5N02g0Gk05dAiPhovuPZ1Cfxdl0d9HWfT3cQq/fBdNMh+FRqPRaHyH7lFoNBqNplK0UGg0Go2mUrRQNCCUUu2VUluVUr87giHeXt821TdKKYNS6iel1Jb6tqW+UUo1V0q9q5T6w/EbGVbfNtUnp3vQUKXUa0qpDKXUr6WOeRVwtbpooWhYlAB3iciZwFDgFqXUWfVsU31zO/B7fRvRQHgO+FhE/gb05TT+XnTQUABeBy4qd6zKgKs1QQtFA0JEUkXkR8d2PvYbQbv6tar+UEolAOOAV+vblvpGKRUNnA+sABARi4jk1KtR9c9pHTRURLYBJ8sd9jbgarXQQtFAUUp1AvoDO+vZlPrkWeBuwFbPdjQEugCZwErHUNyrSqmI+jaqvqhu0NDTiCoDrtYELRQNEKVUJPA/YL6I5NW3PfWBUmo8kCEiu+vblgZCIDAAWCYi/YFCfDSs0BipbtBQTe3QQtHAUEoFYReJNSKyob7tqUfOBSYopZKAdcAFSqkqA0I2YZKBZEdoG7CHtxlQj/bUN66goSJSDDiDhp7upDsCrVJJwNVqo4WiAaGUUtjHoH8XkaX1bU99IiL3iUiCiHTCPkn5pYictk+MIpIGHFNK9XAcGg3sq0eT6htX0FDH/81oTuPJ/VI4A66C54Cr1aa+osdq3HMuMA34RSm1x3HsfhH5sP5M0jQgbgXWKKWCgUPArHq2p94QkZ1KKWfQ0BLgJ06zUB5KqbeAkUCsUioZe1y8J4C3lVJzsIvpVT5pS4fw0Gg0Gk1l6KEnjUaj0VSKFgqNRqPRVIoWCo1Go9FUihYKjUaj0VSKFgqNRqPRVIoWCk2doZSyKqX2OKJ9vq+Ual6DOgYqpZ73cC5JKRVbQ9seUkotcHP8RqXU9JrUeTqg7HzpiEVVWbnPfRXJVFP3aKHQ1CVGEenniPZ5EriluhWIyC4Ruc33pnls7yUReaOu2vMHjqB5/uIS4GcvQs2sBm72ox0aP6KFQlNfJOKIjKuU6qqU+lgptVsptV0p9TfH8ascvY+flVLbHMdGOnNTKKVaKqU+dQTJexlQjuOdysXoX6CUesixfb1S6gdHnf9TSoVXZmTpnoZS6iul1JNKqe+VUgeUUuc5jhuUUkuUUr8opfYqpW51HB/tsO0XR+6AEMfxJKXU40qpRKXULqXUAKXUJ0qpv5RSN5Zqe6HD1r1KqYc92DfHYctXSqlXlFIvOI6/rpRaqpTaCjyplOqnlNrhqOs959O947qBju1YR8gUlFIzlVKbHH+X/UqpRR6+oqmUWv2rlLrW8f3sUUq9rJQyOE5tBiZX9l1rGi5aKDR1juPmMRr7zQPsK2pvFZGzgQXAfx3HHwT+LiJ9gQluqloEfOMIkrcZ6OBF8xtEZJCjzt+BOdU0P1BEBgPzHe0DzMUenK6/iPTBvno6FHu+gKtFpDf2KAg3larnmIgMA7Y7yl2JPQfJIwBKqbFAN2Aw0A84Wyl1fmlDlFJtgQcc110I/K2crd2BMSJyF/AGcI/Dvl9K2V4Zg7ELQT/gKqeglONcYLfDnjOBq4FzRaQfYHVcj4hkAyFKqZZetKtpYGih0NQlYY7QJFlAC+AzZY+Uew7wjuPcy0AbR/lvgdeVUtdjT0xTnvOBNwFE5AMg2wsbejl6Lb9gv4n1rOZncAZq3A10cmyPAV4SkRKHLSeBHtiD1h1wlFnlsNeJUyR/AXaKSL6IZAImx9zNWMfrJ+xhKv6GXThKMxj4WkROOgLjvVPu/DsiYlVKNQOai8jXHmzxxGcikiUiRsfnHu6mTAtH7hSwi//ZwA+Ov+Vo7OHRnWRgj/SqaWToWE+ausQoIv0cN64t2OcoXgdyHE+gZRCRG5VSQ7AnL9qjlKpQBnAXg6aEsg9BpVNkvg5MFJGflVIzscfKqQ5mx7uVU/8/yo0dyst6bKW2nfuBjuv/LSIvV1JHVW0UVnEeyn5X5VOJlv9Mbr9rpVSAiNgc9qwSkfs8tBUKGL2wSdPA0D0KTZ0jIrnY01guwH7jOKyUugpcXjR9HdtdRWSniDwInADal6tqG46hDaXUxYDTqyYdiHPMYYQA40tdEwWkKns496k++kifAjc6J42VUi2AP4BOSqkzHGWmAV97uN4dnwCzHT0ulFLtlFLlk9B8D4xQSsU42p7kriLH953tnFMpZ0sS9l4A2Ie/SnOhsudgDsOeKe1bN9Xv51Sv4QvgSqedjms7OrYV0NrRnqaRoYVCUy+IyE/Az9hDiE8F5iilfgZ+w56QBmCxYyL4V+yi8HO5ah4GzldK/Yh9mOaoo+5i7GP9O7H3XP4odc0DjuOflTteG151tL3X8RmmiIgJe3TXdxzDXDbgJW8rdGRrWwskOq5/F7vIlS5zHHgc++f5HHvY8VwPVc7A/n3uxT7n8Ijj+BLgJqXUd0B51+JvsHsr7QH+JyK73NT7AY5emYjsA/4JfOpo5zNODSOeDexwDs9pGhc6eqxG04hRSkWKSIGjR/Ee8JqIvOeDemcCA0VkXhXl2gBviMiFVZR7DtgsIl/U1jZN3aN7FBpN4+Yhx8Txr8BhYGNdNu7Iy/yKqmLBHfCrFonGi+5RaDQajaZSdI9Co9FoNJWihUKj0Wg0laKFQqPRaDSVooVCo9FoNJWihUKj0Wg0lfL/1+Ewv/c3uF8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(table.index,table['mean'],'-o',color='b',label='Data (no controls)')\n",
    "ax.fill_between(table.index, table['low'], table['up'], color='grey', alpha=.25)\n",
    "ax.plot(table.index,table['sim'],'-D',color='k',label='Model')\n",
    "ax.plot(table.index,table['mean_h'],'-D',color='r',linestyle=':',label='Data (controls for health)')\n",
    "plt.xlabel('Residual income group (e)')\n",
    "plt.ylabel('% deviation in health expenditures')\n",
    "plt.legend()\n",
    "plt.savefig('../figures/fig_d2_spending_by_income.eps', bbox_inches='tight',dpi=1200) \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
